{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "9932a712-28c7-4e05-a39c-ee14239cabc0",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import multivariate_normal as MG\n",
    "from tqdm import tqdm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "db417cc7-7587-4465-a5d4-d793d6ca6a60",
   "metadata": {},
   "outputs": [],
   "source": [
    "def Initial_miu(X, k):\n",
    "    miu = [np.random.randint(X.shape[1])]\n",
    "    while len(miu) < k:\n",
    "        dis = np.ones(X.shape[1]) * np.infty\n",
    "        for i in range(X.shape[1]):\n",
    "            if i in miu:\n",
    "                continue;\n",
    "            for j in miu:\n",
    "                dis[i] = min(dis[i], np.sum((X[:, i] - X[:, j]) ** 2))\n",
    "        miu.append(np.argmax(dis))\n",
    "    for i in range(k):\n",
    "        miu[i] = X[:, [i]].copy()\n",
    "    return np.array(miu)\n",
    "\n",
    "def KE(X, miu):\n",
    "    cluster = np.zeros((0, X.shape[1]))\n",
    "    for i in range(miu.shape[0]):\n",
    "        cluster = np.vstack([cluster, np.sum((X - miu[i]) ** 2, axis=0)])\n",
    "    return np.argmin(cluster, axis=0)\n",
    "\n",
    "def K_means(X, k, echo=False):\n",
    "    threshold = 1e-6\n",
    "    n = 0\n",
    "    \n",
    "    #初始化均值\n",
    "    miu = Initial_miu(X, k)\n",
    "    \n",
    "    while 1:\n",
    "        temp = miu.copy()\n",
    "        #计算类别\n",
    "        cluster = KE(X, miu)\n",
    "        \n",
    "        #更新均值\n",
    "        for i in range(miu.shape[0]):\n",
    "            miu[i] = np.average([X[:, [j]] for j in range(X.shape[1]) if cluster[j] == i], axis=0)\n",
    "        \n",
    "        if np.sum((temp - miu) ** 2) < threshold:\n",
    "            break\n",
    "        \n",
    "        if(echo and n > echo):\n",
    "            break\n",
    "        \n",
    "        n += 1\n",
    "    \n",
    "    return miu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "5d77bd2e-bc8a-48b3-b71d-77e58f1175df",
   "metadata": {},
   "outputs": [],
   "source": [
    "def GMME(X, k, alpha, miu, sigma):\n",
    "    gama = np.zeros((X.shape[1], k))\n",
    "    for j in range(X.shape[1]):\n",
    "        for i in range(k):\n",
    "            gama[j][i] = alpha[i] * MG.pdf(X[:, j], miu[i].flatten(), sigma[i])\n",
    "    gama /= np.sum(gama, axis=1).reshape((X.shape[1], 1))\n",
    "    return gama\n",
    "\n",
    "def GMM(X, k, echo=100000, echo2=False):\n",
    "    #X(m, n)\n",
    "    #初始化\n",
    "    miu = K_means(X, k, echo2)\n",
    "    sigma = np.array([np.eye(X.shape[0]) for _ in range(k)])\n",
    "    alpha = np.ones(k) / k\n",
    "    \n",
    "    for t in tqdm(range(echo)):\n",
    "        #E步\n",
    "        gama = GMME(X, k, alpha, miu, sigma)\n",
    "        \n",
    "        #M步\n",
    "        alpha = np.sum(gama, axis=0)\n",
    "        sigma = np.zeros_like(sigma)\n",
    "        miu = np.zeros_like(miu)\n",
    "        for i in range(k):\n",
    "            for j in range(X.shape[1]):\n",
    "                miu[i] += gama[j][i] * X[:, [j]]\n",
    "            miu[i] /= alpha[i]\n",
    "            for j in range(X.shape[1]):\n",
    "                temp = X[:, [j]] - miu[i]\n",
    "                sigma[i] += gama[j][i] * temp.dot(temp.T)\n",
    "            sigma[i] /= alpha[i]\n",
    "        alpha /= X.shape[1]\n",
    "    \n",
    "    return alpha, miu, sigma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "4f8c3d9d-97cb-415c-bda1-1da092abef67",
   "metadata": {},
   "outputs": [],
   "source": [
    "def draw(X, Y):\n",
    "    color = ['r', 'y', 'b']\n",
    "    shape = ['o', 'x', '+']\n",
    "    for i in range(X.shape[1]):\n",
    "        plt.scatter(X[0, i], X[1, i], marker=shape[Y[i]], color=color[Y[i]])\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "a9082003-d01a-415a-bf7d-d6f0cbdfceae",
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 300\n",
    "num = 1\n",
    "miu_real = [np.array([-num, -num]), np.array([num, -num]), np.array([0, num])]\n",
    "sigma_real = [np.eye(2) * np.random.rand() for _ in range(3)]\n",
    "alpha_real = [0.2, 0.5]\n",
    "Y_real = (np.random.rand(n) > alpha_real[0]).astype('int') + (np.random.rand(n) > alpha_real[1]).astype('int')\n",
    "X = []\n",
    "for i in range(n):\n",
    "    idx = Y_real[i]\n",
    "    X.append(np.random.multivariate_normal(miu_real[idx], sigma_real[idx]))\n",
    "X = np.array(X).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "c436d002-b338-40c8-bbc6-96ad0ae8030e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAArV0lEQVR4nO2de2wc15Xmv9ts2hbftuiR2HQovuQFgsCiJUayDVvUYrI7nkViY4INkoGCnZn8oQ2wFiVLi8FmhLUkzBrYQRK9nJnFChtjHitMZoMkG+9skEkCxJKz0SOiRGeSOLFEUo7dLT8kSxRJySa7++4fxVusqq7qru6qrld/P0Cguroe91Z3f/fUueecK6SUIIQQEl9SYTeAEEKINyjkhBAScyjkhBAScyjkhBAScyjkhBASc9JhXLS7u1v29/eHcWlCCIktExMT16SU91u3hyLk/f39OH/+fBiXJoSQ2CKEeMNuO10rhBAScyjkhBAScyjkhBAScyjkhBAScyjkhBAScyjkhJDEsW2b9q9RoJATQkjMCSWOnBBC6oGywk+eNL9++eUQGhMgtMgJISTm0CInhCQGZXk3iiWuoEVOCCExhxY5ISRxNIolrqBFTgghMYdCTgghMYdCTgghMYdCTgghMYdCTgghMYdCTkhINFo9EFI/KOSEEBJzGEdOSMA0aj0QUj9okRNCSMyhRU5IwDRqPRBSP2iRE0JIzKFFTkhINIIlzqeOYKBFTgghMYcWOSHEdxiZEyy0yAkhJObQIieEeMLO2mZkTrDQIickQTDtvzGhRU4I0anGgnbjB6clHgwUckISACcXGxsKOSGkpoGAfvDoQCEnJAFQVBsbCjkhpKqBwLoPB43woZATYiHOVm0c20y8QyEnhOi4scQ5oRo9PAu5EOIjAP4WwBoAEsBxKeVRr+cl3uCPrHooVCSu+GGR5wHslVJeEEK0A5gQQvxQSvkrH85NSKIJYrDweg3r8RzgoodnIZdSXgVwdfn/c0KI1wD0AqCQhwCtytqhUJG44quPXAjRD+BhAGdt3tsBYAcA9PX1+XnZ2ELBaFyCGHC9XoNGQXzwTciFEG0AvgVgt5TylvV9KeVxAMcBYHR0VPp1XWKGVqV3eM9I3PBFyIUQzdBE/ISU8tt+nDPJ0NIhQQy4bq/h9D6NgvjgR9SKAPB1AK9JKQ95bxLxA/7oCGkchJTevBxCiMcBvALgnwEUlzf/mZTye07HjI6OyvPnz3u6bhKgpRMevPelT4ZjY9rfRr4nUUcIMSGlHLVu9yNq5ScAhNfzxAUKACEkajCzM0Q4GARPNfMTSR+06QNPDhRyl3CCMt74tWpOFD73KLQh7iTtHlLISUPhxgpttEHbr34l/T5FGQq5S/gYGk+cRLlaJie1Y8MQd2ubozDAxPV3kNRBmkJOGpJqV75JyoLGUcwgJd6hkFcJv5zxwq8nqTCeyJzCA8MME7QT7clJYGQk+LbUQlKfrCnkJPHUKjZJ+ZED9bWaR0bMTy1Jum9xgUJOEoWTmCix8UqQIhVF6zGKbaqFuLbbCQo5SSz03a4QZG0XEjwUchI77MQo7r7bclgFMgoDEkU7WlDISWKJku82Cm0I6/pR6XuSoZCT2FDOVZIU32056CoiTlDISeKJgiXeiOLbyH0PGgo5iQ1urO6kioRdcpLXhSPiQJzbHiQUcuIL/MHZ0wguHyfq1fdGvJeVoJCT2BH3H7BbITJG3thF41RyXcTZtRHntocBhZx4gj84d3i5H5OTfrUiHPz6LoRZuCzqUMhJ5EjqD9TtoGfdr7MTaGoC2tq017Oz2nuVaq84uTbicH+TXLisHlDIiSfiJA5xZXZW+zs/H3/r3Av8rjlDISeRIeluGrdCpLZ3dWl/lZA//rj9ftVeN073N8ptixIUcuILjfCDC1r4rOUFlHshKaUHaqURvmvVQiEndaNa4Uvao7NTP2q9H4paKzl6ub9J+UySCoWcxIowBCVsl4SyxBmxQZygkEeYuP5YvQpfUP2t1/0NW/grYRwY3MSyA9HtC9GgkJNYEKagRMHlE4U2EO/U6/OjkEeQuFtBbkTHzz5Ve656319j/1W4YFQ+u2r7zgEkHlDISSyIgqBEQcSi0IZ6k8RBo97GA4U8gkRBtPygnCXuxxe61nMFcX/VuVUWZlQ+y1r7Hna7SXko5CRWUFCSS9xdiuWot/FAIY8wSfgCW/H6hTYe5/Vc9by/QT9V0cJubCjkJPEkybLzSpTvRVJciuWoV58o5CQUarXEy63XGTRuBCcoSzyJ7gjiHl+EXAjxIoBPAnhXSvkxP85JiFcocivE6V5EsU1Rxy+L/K8BfA3A3/p0PkJMROmxO0qiGKX7QsLDFyGXUp4SQvT7cS7SmNRDiMISOT9CKv1ua1IFP8j+RPne0UdOYoXxRxTWDyuKohiFNpDwCEzIhRA7AOwAgL6+vqAuSyJOEG6KWs9V79R/u2XM6u2uSYrgB+neipIrzYnAhFxKeRzAcQAYHR2VQV2XJA/ryjlhW+aEhA1dK8SRIAQyim6Keqf+lzu/GqSicB+iTLX32sv9jOJ31Ipf4Yd/D2AbgG4hxFsA9kspv+7HuQlRGOuXANrq8kA0f1hRJsqCRGrDr6iVP/TjPCQahOETjJKoOC2xVolaLXdAK3e7bVv47qK4UctTj9/XigJ0rRDPBCU6cXjEDQKvdWqCnCBs9LDAoKCQNyBeFhNoxB9NJeHzKpDqOEBzF9W6uLIfJOXzbbRBn0JOaiassKyk/yid8Ot+j43Vdpwb6v2dCCNkMw5QyBuIar/4dpa48djJSc16TCpurbpaIyiiZDUmVRTj3n63UMhJzSgXQJR/9Na2Rbmtlai1D1aRVhZ5PajXfS430MT5M/ULCnkD4eVHlgQhrBWvizVUsnbrGQPtFqfPtxE/7zhCISee8ftH7od4WMUzKtmgflBtm8MQ5SCLfsXxM/QbCnkDktRY2nLMzwNtbeFcu1ohrcVf7ZdIW9uaNJ95UomVkPPLlGz8FA+jIKlJ2Ub3pyahz0noQz2IlZATUg1KxGdntcEh7Eiban3t1VjiflvO9JHHi1gIOR/z/KOecb1eqYd4jIyYE27CTLYB+N0l9SEWQt7o8MdfG3G0Ko1tddNeN330w0VFok0shDyOP8io4fdTTT2fkpL4+dpF0YT9dECSQyyEvFGhS8kf4nC/vH7W5Sxxfn+ST6yEnF/A2vH7qSbOT0lhtFlFzHR2apOvxnh2P9oRx88hKBrh3sRKyBuNOItlFFH3MYrU47Pm96dxoJA3GPXKuIsDk5Pa37CyO+shrHSfONNI94ZCHgOS+MULEvUDVgIeZerxWfP7k3wo5KThCHutTz+vy2JXzjTSvaCQk8Rj/UGT2qi1fG6SBTQqUMhJw5FEYWGxK2caoe8UctIQUNBqp9rBgYNJ8FDICUkQUfALq+ggEhwUcpJoaB16p9rBwVpdkve6/lDIScMxOelfRmVUCaNv1kFTRQeR+kMhJ4nGzppk9EptVDs4hFX3vRGhkJOGQVnidLPUhyj45xuVVNgNICQIXn6ZFiJJLrTISSzww8qjxRgMvK/BQ4ucENKQbNuWnPkSWuQk0tQjfLCeFiOtfRIGFHJCaoCCHV+SmFtAISeRJi5+7SiLQ5TaQuqDL0IuhHgSwFEATQD+h5Tyv/pxXkKihpNgk/gQF+OgGjwLuRCiCcBfAvhXAN4C8DMhxEtSyl95PTchiqj/2KIoDlF+SiD+4odFvhnAZSnlNAAIIb4B4GkAFHKSOKIo2KQ2kvTZ+SHkvQDeNLx+C8AW605CiB0AdgBAX1+fD5clJHpESRw46DQOgU12SimPAzgOAKOjozKo6xJSDyiKJEr4IeRZAB8xvH5geRshJAJw0Ek+fmR2/gzAeiHEgBDiLgCfA/CSD+clhBDiAs8WuZQyL4R4BsA/QQs/fFFK+UvPLSOERB4pJYQQjq/pnw8GX3zkUsrvAfieH+cihMSDmZkDyOdvYnj4MIQQkFLi8uVnkU53YWDgQNjNayiY2UkaAlqG3rBa2sViEfn8TWSzRwEAw8OHcfnys8hmj6K3dxe2bZMABGPYA4JCvgy/aITYY2d5T03tQVNTJ3p7dyGbPaoLem/vLgwPHwYgyp+U+ArL2FbDiRNAfz+QSml/T5wIu0XOxKmtdUSVKj15UvsXpdKlUWqLE1JK3fK+fPlZ3X2SzR5FoTCLoaFDpv2V2L/8smYUjY1p/9RrUh8a3iJ3ncZ84gSwYwdw+7b2+o03tNcAsH17fRtZLXFqK4k0QohlCxsllvfQ0CFMTe0x7X/58rO6mJPgEFIGn5szOjoqz58/H/h17bAK+diY9rdEyPv7NUG0sm4dcOVKXdpWM3Fqa0BEyXXm+jsXIaSUOHly5QF+69YCpqb26D5xq4+cYl4fhBATUspR6/aGt8hdpzH/9rfVbQ+TOLU1ZKIk8HYYJxm1tkq8/HKwAqncKUampvYgne40ibay3NPpLtciXil8kbij4YXcNX199lZuFOvGxKmtAREloXZrPFgnGQGJO3emMDPzPwML7zP6xO0s76GhQ7rwKjF3K8QMX/QPCvkyFX/ozz9v9jsDQEuLtj1qxKmtIRH1Eq/GScbPf/4ZrFo1hJMnBYBhfO5zn8eqVcFY5kIIpNNdjpZ3KpUq2d8Nxv4BpeGLflnmjWLxU8jdoiYJ9+3TXBR9fZowRnHyME5tbWDKDRpGwfzww7fw4YdvAdgGAFi1aghBhvcNDBwwCWC1lrcd5SZR/fKvN5LF3/CTnaSxCdUSP3Gi4mBrnGTcvfvH6OoaC9xHXg6vFq91EnVsrKiLrhcxr+QSMg4WcbLaOdlJSJRwESJqN8l4584UpByKhNB4tXjt+nf58rN6WKMXy9mtxZ8Uq50JQaShCS1RZd8+8xwGoL3etw9AqUU5NlbEN7/5v/EXf7FeT8yxPk0H+XRdLlEon79ZsS3G/VtbRwAAra0jyGaPYmJi0/J5bujnsetvJYxirjAKtrEPly7trroPUYIWOYkeyuXwxhtAUxNQKGgx8Eny81cIEa00yXjlysFQLUkhBIaGDkFKabJ4M5lxVz5uIQTm5yfR2jqCjRvPY3p6r36OhYVJNDevNYn4pUu7MTd3Fvfd96Tr/jlZ/Ol0J/L5WQwPH8bw8GFIKZHLHUMudwyAv376oKCQk2hhdTkUCtrfpGWnuggRdZpkBKBbjkD9Ij7KMTNzAO+//310dJhXdZybO4srVw5WFFspJdraNAt8enovhoYO6f0BgKWlt5HLvbA8SEEX2Y6OLa76V85H3to6goWFSQDKQjcfGzcRBzjZSaKGU1aqIinZqdYBC9BCRI8fdzVQGYVKEdQknrKQlbha6e0dx/DwkarE1opRbKs9r8LJ/93U1IlCYdb2utp1omuRO0120kdOKhNkAa5K2adJyU7dvl0T7XXrACG0vy5FHCjv/wU0EVO+a2BFNGdmDtier5y/3c7YGx4+jJaWDQ7ncuevV+4ZI1u3FtDbu6tExAFgcPCQLshO7TIyMHDAdE/UPRscPFhy7zKZcYyNFfVqjsZ7FwfoWiHlCboAl5PLwfh+Uti+veZ76OT/VQJVTbLNzMx+3WcshECxWNSjRtS57Hzx3d1P47e/fdXUhkxmHOl0l358f/9+x6cCVQ7XyNTUHgwOfhXvvvsNLC29Y3rv9Om1WLPm3wOY161quzkBqzvK7rX13qkm1lJmIArQIiflqRBd4TvPP6+5GOxgdioA+4gWoyUJaIKktp08mXIsZjU9vR/Xrr2kH1ssFk1RI0tLNxwiU7T37FCDyPvvfx+XL+9GsVg0tXt6en9JH5Qlns0exYULm3QRF+Iew3mvYWbmK8hmj+L69Zdso0sqPYnY37txZLPHTPcuTqGHAC1yUomgC3AZs1KTHLXigfIRLZ26UA8PHzb5ga0iLqVEoTCLhYVJPfRP7d/aOoKhIW1/IYQlFnscxaLE1asv6K9v3DiJ27df1f3mPT07kUoJZLPHcPPmKWzaNKFXS2xtHYEQcOzD/Pwk2ts3486daeTz10x9v+uuDwFokS3WWi9SSiwt3dDb4PQkUnrdIwBE7KxwI5zsjCsusgJ9gSVxI4vVRTI9vR+FwqwhsmU3stmVCUk7i9xpwnHr1oJeR8Uu+/LKlYN61Mrw8BFIKXHqVJNpHyklJiY2mfzdahLTKMJ2rhcpJYrFIl55xdnW7Ot7Tu/vistkN27dOou5ubMV+x2XbE4jnOxMEspv/cYb2syS8lvXYxLSztVBF0cksLOuVXKLUcQzmfFl90HpJJ7dhCOg+apV0oydL76/fz82bjy9bM3CdoEJIQQ2bZowbTeK+NTUHly5csCxf9brWlHuFbPb51hJSGQtEShhJlvVAl0rcaSc39pvq5wFuGKBNSVdkcmMY/36I2qvEveB8okbUW4WQOox3NZYbAAGy98+XluzckvbOjj4VZObZd26/UilUrpFnU534erVr2NxMYumptUoFK7rx+bzTUinC7plb3UJ9faOw6q51lWLKqXlxzFtnxZ5HAnDb33lClAsan/DEnGuQ1oWu5DE9euP6O4L6ySeihpRgmgM/WttHUE63YXm5ntL/Ni9vbv0AcHJX5/JjGNu7iyy2WN6Cr7ilVfSyGaPorl5LRYWJjE19awu4tnsMVy//n1ImQcAFArX0dKyAc3NawAA6XQBLS0bsHr1p9Dbuwvd3U+Zzm0ceOzCCbW0/JUJ3GKxuHxdbeK0UCh4Kj0QFvSRxxEnv/Xq1cC1a6Xbk4DHBJpGwE2SkBVlfQ4NHTJYxVrSzODgQf28lfzJdvvMzBzA9esv6e6UwcGvmnze6XQP8vmrJW1qbs7gkUd+izNnHsDS0tv69paWDejs3Irm5nv1rFdl2Sva2jajo+MRDA8fRiqVMoVSKmt7aekGhIBp/qC9fQsefvin+spH+fxsVfcxKJx85BTyOHLiBPAnfwIsLZm333UX8OKLyRQ2TrqWpZqyrXbH1mviT8WoK5+4URzT6e6SqBRAc+0o37pxAnXr1gKEEMt1Zm7o1ndb22ZIuYilpXexuJhDJrMTyo2kxZp3YmDgoOke9fTs1KNuAHOUjfLhWydvwxZxgJOdyWL7dqCjo3T74mL94ru94IdLhOuQlsXJxWF0g5Q7ttxrayal28xKABgYOGgScWOlQ03Em0z7p9Pd2LjxvB7xYmRiYhOKxeKy6+MY5ubOIpPZCSkXsbAwiWJxEQDw3nvfRC53zBBrPqsPTsrtYxRxALh69QWTiNtN3kbVrQJwsjO+vP++/faoCJuxgqEQ0Gegas0M5TqkFanHSj7K9aIyKZXINTV1oFC4pbssyln1qVQKTU3aQs0DA1/BzMx/NFjmBdP18vlruHjxUXzwwQzy+Wt6dcTTpzV/+oULo9i4UXuaz2aP6mGGxtosyh2jXDrWe1DudhgnYp0meKNgmVuhRR5XnAQsCsJmDI8EVkRcUUtmKMMgXVHJuq4GY71uZd2qrM9c7r/pE4DahOFK9qTKrjRmdBYKs3jvve/g4sWPY3Dwq2Wv+8EHbyCfv6Zb59PTe/XXTU0daGpqKpnUtYY5KuyE99ats7b7AsD09F590Kn2ySZMKORxJcrCZhceaaXaJwePRaZIbSgRU9au+qusZaPbRIm6En/lClEW7dLSu1hYmMTp02stVxHo6XnGVISrtXUD8vlreoRLb+84Hn30bWzceBKFQgGXL+82ncHqhlHYpevPzZ1FT89O9PaO6/up19nsUf3Jw/pkE9XQQ4CTnfEmqOzOakmlSq1wK5ykjDTGWGoApsxOO4wWrIpNt8voFOIeSPmB4cgUAM1yb2vbjGLxQ+Tz17C4mDWdv6dnJx588ChmZg4gl/urZQt9DR555C2cOdOz7G+/G8CHphhzq3vF2C/jpGlz873o798f+XhxRq2Q4KhUU5xhg5HGGN2RyYyXhOrZYY3qKBaLpqgPAMvnakc2u/LU+NhjH2JycjNu3zZXUbSrR7527TOYm/uJaXtLy4blY7UBobl5DbZseQsXLowine5Ee/vDJmFWeqfaWiwWTWUCop6qz6gVEhx2bh/146BLpO54TS83RnfkcsdMIq6iTtLpbtMxdi4Mu3bduPGPpm0//enduH371ZLa5polPY6enp1obtZcMW+//bUScV8ZADSr/v77P4upqT24fftVSPkhhoYO6SKufPfG9mhlAg6a+h5HPAm5EOIzQohfCiGKQoiSUYLEDL8yJ+382X/3d5q7JczM0AbA7YISTmJvtFhXUvs1envHsXHjeT10sLl5Lfr6njNlTxYKBVy48AhyuWMl4nz16gsllrfCziJXFRa7uz/juv+53DE9tLCjY4vJ0q5HxmZUarJ4DT/8BYBPA/jvPrSFhInfC0h4WDSB1IZRrADnMq5OtUTm5yfR1jaiT/RZJxRv3Di57Ia4CwCwevW/RX//fv39ubmLuHjxMXzwwRUA0C3txcW3kc+bF4lYcYmskMlok4+53DEsLExiYWFST9Splt7ecb0ML1Bai2alNkvtGZtRqsniySKXUr4mpfyNX40hIRL0AhLEd4yhck4LSpS3TGf1aJNLl3bpLhVV6+T27Vdx6lQT5ufPoa1tMxYWzuHChUcxNfUsBga+gkJhFvPz55DPvwshWgFoYp7Pv4OmJrMrxskyt4YVqmxLAGht/XhJ3RYn7BaUqLQ8XjXUy8KvFfrIiUaSMycbqNhWJbFSZWvtxF65TRYWJpHLae6J1tYR3L79Kn7ndz5bcp25uXOYnz+HbPYYLl78OBYWVsRZygXT/oXCNfT2jqOv7zmk02tK2p1Or0E63VniW5+dPYNUShsUOjoewX33fQqp1H36+z09O9HUtLrkfLncC8jljmFp6Qamp/cbCmSVluRV8e4rbXe33qjbVZiCoKKQCyF+JIT4hc2/p6u5kBBihxDivBDi/HvvvVd7i0l9iHKCkRfsard/4QtAd3cihd2pfrgSp5mZA5ia2lNSg7ypqRNNTU0O9cNLS8PeuTNTsl8lisUilpZu6G6W5uY16OnZCQDI59/B9ev/R6+doizv+flzKBa1QeH69W+jUJjF2rUrLrurV1/Qy9xarf7m5rUYGjqk12lXyUzGMgHGeHeg8iLVRvy08L1SUcillJ+QUn7M5t93q7mQlPK4lHJUSjl6//33195iUh+inGDkBTuX0eIicP16/RflCJhKa3kak3XOn99oOvbate+iUCiU1BhR583ljumlbjOZceTz77pul5r0vHr1a6YaJ0tL7+Dq1Rf0hS/S6U5kMuPo7HzEdmBYXMwilzuGXO4F3Z+u6Ol5Bo899rZpWyazQ88CtSY1WZ9A1EIa1bhHKg2aQULXCtEIOnMyKHeHG9dQQuYCKhXOSqVSGBo6pPu7FWri8fTpjMliVbz77v9CS8sGXfCGhw9j1aqHyralpWXDchVC2IYXGhkePozh4SN4+OGXsX79EQwOHtJDDo3nM6LW5VTcuvUKLl58zLTNWizLek31BFKLe6TSoBkrH7kQ4g+EEG8BeBTA/xVC/JM/zSKh4McCEm4EutxSdX4LvFvXUBLmAqAVzrL6xI3p5UII3HNPv+mYe+8dAwAUiwtIp9foWZFPPJFfDjXUXCHt7VuQzR7FqVNNuHPn52Xb0dU1hqGhw2htHUEq1YbFxbcd952Y2GQKf7xwYdRUhxxwnhxVbVxYeBXz8+fQ07PT9knEyXJOpVI1uUe8VJusB16jVr4jpXxASnm3lHKNlPL3/GoYCRg/BNTtWqJOETK7dvm/Fqmdy8iOGM4FOMUwVyqc1da22fRareJTLC5g1ap1uvU9Pb1Xdz/cvv0q2tvNa2E60dy8Bul0F6an92JhYRKp1D3I599Ba+uI7hMHNCt71aqHsLAwiYmJTbprZ2FhsiThSJ1XWfmKV15JY2FhEs3Na9HWthkPPnjUIqqdpmqG1Yh8JSoNmkHCFP0kUm0NFr9W33G7+IObWizljldtdtvHEyeAP/ojoFCwfz+GJQNqiWEut7ACoMVeDw8fQbFYxPT03pK1P+1S9TOZnRgaOozp6T3IZo8hlWpBsbjyPdIs8nbcufOaaRGJ5uY1uOeeddiw4f8tR7xMmo5ZWJjUM0sVygVkjDdXPPFEHqlUyraUrtO9UqV5a1mMw4l6LtIBMEU/2Rit6e5uLSqjGqvWyUL+/Oers87dhjBWa/1aj3dr+Su2bwf+5m/sLfPVq2Mn4rXGMCt3QCYzbptkMzSkuRimp/cine60HKuJeHPzWmQyO3HXXRkAWpLQxYuPQYgOpNPdJhEHsLzgw1zJSkD33/9ZzM2dw9mzH8F9933K9N7q1Z/SBw4j+fy15cHG7ApRbbbeIyWgAwMHbKsZDg4e9NU94jarth7QIo87dta0HeWqDVaykN1arG4tcqcngFWrtGiSSsfXuuxbVKtF1kAt63OuHLdbd6dYLeGOjidw9eoLJe+1t29GofCh7qtubl6DpaWVbE1lLVstcoV1WTe7zE6F0fpfsZS1NtstyWZdCEKrauj+acUPK9rLUnvVwOqHSaVSpUGFENokZq3ncFN2thoXjZ2oAu6Odxp4yvUxgUgpTeVl3a4rOT2937Ao8rhe00ShRLalZQO6u5/W3Q/asmooccm4pbV1gylpyMrWrQVdlNvbt6CjYwuGh4/opXHVIso3b76MfH4WmzZN6AssT0xsQjrdiZGRH1cUVAB1cX/UOrhWg5OQc6m3uOM22qKcO+P55ytb9W6uo8TWjdVbrhZLpeO57JtjDLMb0RgcPAghVGTJIUhpFvLVq58CoEWK3LhxD7q6/jVaW0eWa3YfcBTytraPQ8ol2xjw1tYRtLc/XlbIjYlK6XQX+vv3m3zcKzVgbixnkz6GjRtP4/LlZ/XYcAD6OexqqszMHEChMFuX+ijKNWMU8qAShCjkccdJ1IxUSuxRQrlrl71rQ13HDV6LZbk53m7gSULykkvKPcYD7sRjYOAgpqf34+LFx1AofGB67803/1z/f0fHFhQKs1hYmERn51ZcurSrTLuWTAs6GN0zqgiWFbUmp3FydWjoEFKplN5XayEw9TA2N3dWfyLJZFZ855rlbvbxa8dJ/UlEbbMrKlYrXgZXr1DI446dqDU3Ax0d2gLNbn3B27drlrCdkAsRLZGsxvJPIE4xzABcTdIpd2qhcFNfvLilZQOEgMlizmR2Ynj4iP7aGCXS1LRaT41XKPFevfop3dqfmnq2JNKlt3ccs7NnMD9/bjnMcY8+0aoSl4x9tataaI1oUV1WwmxNarp0afdy/yb11Hw/KiAq/BhcvUAfeRLwaxKv3KRnSHWWiTPKgnT6a4d1CbcLFx7VxdxKJjOO9euP6Oc0+uPb27dgbu4sWlo2IJW6e7mI1lm0tW3Gpk1n9P0uXdptEtzW1hFs2jQBIQRef30X5ufP4b77nsTAwIGy7bZeX1tf037VImP4ojVkUpW3Na5e5HZuoRJBlLWljzzJ+FX728lNs26d93MT3ykXI20nHvZuCucBWgnw8PDhkvR3KSVaWzdg9eqn0N9/QK9fnk7fa2pHLrcSeaJEXaX5P/jgUb0fxr/q/MZFIS5d2m26vha1osW9WwcL9USwYsmvvDc0dLiknoxf7g/rYKSeJugjJ8HS4L7nuOF2IQmFk5vCDmXVzs2dxaVLu3Wrvb19M4rFRczPnwMAdHZuxdTUs0in79UjTNS1rO4fZd1Xcv9YJzeVUCtrXz1FSKkE33x8oTCrT3hafdZqUeh6uT8qZdXWCwo5WaHBfc9xo5ZVb+wiK+xob38cXV1bkU7fC0Ct3iP1OuWAluovxEpMt5VaLFTj4HTz5kls3HheH0SkXMT09H5s3Hgaly7tRjrdtRyuuGKdrwizhJQwPREY/edqkKhmbiHK0EdOSMypJp7cLtZ51aqHTEWw2to2Y/Xq39fD/9RxAEzXUajEHOMkpRdUXLhd2n5r6wi6u1fcOeVcSwAc3uvEwMBB0z2Ji4jTR05IAqkm5M0usuL113eVxIV3dj5iEnHjee1Q2ZWVJvXcTs6mUik8/PDP8JOfNOvbjOGMXV1j+vZyVr+xiJj6v/E9O998XKGQExJTqg15s/qtAeg1V9rbt+jJNdrxosSSVcW2bt16xWQtG/3OTtatspzT6U7k87O6+KvCVcZBQGWeWjH6tq39Mt4TaxKRlrJ/A1JiOalpf2iLJNcLCjlxT4JqlSSBWuLJrRZsOn2v7l92Ot5YbMsYi21M9jH6na0Y/d7quJs3T5rOowYBY9KOtT6L1v6vVAytHBo6ZPKzp1J3YW5Om5zNZMb1ui1+JAFFBfrIiTv8KnVLfMdr0Se3x0sp9YJUQ0OHTLHYfX3PYXDwYMkxxmOtvnmF1cqent6PXO6vSkQcMMehW9tsfDoZGjpU4mcvd824wKJZxBu1VhwkiUEJvCpgVUvlRbvJUrvJ2amp/4w33/wv+utMZidu3jyl1yMXQnuaMLpGyg0Wbq4ZB1iPnHjDba1xkkhUrW2jiLe2jqCv7zlXa1XaTcoqrMdJKVEszpn2yeVeQFfXGDKZcczNnUU2e6yk9rrRNVQJN6sAOa3AFEXoIyfuYMXBhsWaeJROd5oiSLTkG2nyq1szM40x3NbCWsbJWQAlLhKVIWrM3rR7ArAbLFpbR9DZ+YQe/76Stl8+CSiIdHs/oZATd9Sa9ckJ0shTyUduTTxSGBdykBLo79+vH28UPW2ytFNfQ7OrawyDg1/F9PReNDV1oKtrTB8EpJS2GaGAuWiXk4irAaCpqVOfNE2l7l5e61PoUStAaYapMRyymozZKEAhJ+6oJevTOkGqlmQzno+EilvL06nWNqAl3eRyx/R9rKKnTZDO6klDyj0zN3cR7e0P69uN1y4NMTS32xorbxfB09+/X49vV4OM2t96fut9GBo6hJs3T/peJbFeUMiJe6otzuW0Fui+fRTyCFCN5Vku8ahcmQAA+jWklFi//ojuY29p2eB4bWMbte2lqfbqONVGa2hlKpUqW6qg3H2YmtpTEvESVREHGLVC6gmXZIs8bpYnc7MeJQDHMgFSypJyucpXfevWWdN2O6s3CH+13X2wrlsaBYucUSskeJwmQjlBGhnsIj0qZYSqY9Rq80Bp+r41KqSjY4vlulp5Wet2O6EcGDhQ4kYZHj7s66Sj3X1QiUpjY0VXkTlhQiEn9eP557UJUSMsixspnFwmVrFyElOV7q6scyfRs2qfqlho3e4klPUuD+sU8WKskqgGrii6VyjkpH5s365lfq5bp5lg69YxEzRCWF0mlSxPOzF1a63ncseWS+GucOPGSb3MbJhWr9N9WFiYxNTUHt3v7vdTgJ9wspPUF79WLyK+U0utFjsq1R031mkxks9f0xdNDrM2uNv7EEVLXMHJTkIaHK+1Wtyc31ioSpss3Q3jghB2iURBU+/74AesR04IsaXe/mfN4r3XYvEegV1STsgRIWVfRxkKOSGk7oS5MHEjwMlOQkggxNnijToUckIIiTkUckIIiTmehFwI8WUhxK+FED8XQnxHCNHlU7sIIT4Sp9rapHq8WuQ/BPAxKeVDAF4H8CXvTSKE+IlaFGIly1JLgJmZORBuw4hveBJyKeUPpJT55ZdnADzgvUmEEL8wVvZTYq6yGK0r7JD44mf44RcA/IPTm0KIHQB2AEAfiyYREgjWRSHiUFubVE/FzE4hxI8ArLV5a5+U8rvL++wDMArg09LFEM/MTkKCxbrwcVwXH250as7slFJ+osKJ/xjAJwH8rhsRJ4QES7lFISjmycBr1MqTAP4UwFNSytuV9ieEBEu1FQ5JPPHqI/8agLsB/HB5ZD8jpfyi51YRQnzBrwqHJNp4EnIp5bBfDSGE1AfWOUk+zOwkpAFgnZNkQyEnhJCYQyEnhJCYQyEnhJCYQyEnhJCYE8qanUKI9wC8sfyyG8C1wBsRHdj/xu1/I/cdYP9r6f86KeX91o2hCLmpAUKct0s5bRTY/8btfyP3HWD//ew/XSuEEBJzKOSEEBJzoiDkx8NuQMiw/41LI/cdYP9963/oPnJCCCHeiIJFTgghxAMUckIIiTmREHIhxJ8LIX4uhJgUQvxACJEJu01BIoT4shDi18v34DtCiK6w2xQUQojPCCF+KYQoCiEaJhRNCPGkEOI3QojLQoj/FHZ7gkQI8aIQ4l0hxC/CbkvQCCE+IoT4sRDiV8vf+11+nDcSQg7gy1LKh6SUIwD+EcBzIbcnaH4I4GNSyocAvA7gSyG3J0h+AeDTAE6F3ZCgEEI0AfhLAL8P4KMA/lAI8dFwWxUofw3gybAbERJ5AHullB8F8AiA/+DHZx8JIZdS3jK8bAXQUDOwUsofSCnzyy/PAHggzPYEiZTyNSnlb8JuR8BsBnBZSjktpVwE8A0AT4fcpsCQUp4C8H7Y7QgDKeVVKeWF5f/PAXgNQK/X83pdIcg3hBDPA/h3AGYB/MuQmxMmXwDwD2E3gtSVXgBvGl6/BWBLSG0hISGE6AfwMICzXs8VmJALIX4EYK3NW/uklN+VUu4DsE8I8SUAzwDYH1TbgqBS/5f32Qft0etEkG2rN276TkgjIYRoA/AtALstHomaCEzIpZSfcLnrCQDfQ8KEvFL/hRB/DOCTAH5XJiy4v4rPvlHIAviI4fUDy9tIAyCEaIYm4ieklN/245yR8JELIdYbXj4N4NdhtSUMhBBPAvhTAE9JKW+H3R5Sd34GYL0QYkAIcReAzwF4KeQ2kQAQ2hp7XwfwmpTykG/njYLxJ4T4FoB/AaAIrbztF6WUDWOhCCEuA7gbwPXlTWeklF8MsUmBIYT4AwAvALgfwE0Ak1LK3wu1UQEghPg3AI4AaALwopTy+XBbFBxCiL8HsA1aGdd3AOyXUn491EYFhBDicQCvAPhnaHoHAH8mpfyep/NGQcgJIYTUTiRcK4QQQmqHQk4IITGHQk4IITGHQk4IITGHQk4IITGHQk4IITGHQk4IITHn/wM/b9SRbJpAkgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "draw(X, Y_real)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "2dc75ba8-87fe-453c-87f6-e7e588f00387",
   "metadata": {},
   "outputs": [],
   "source": [
    "k = 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "56ae933f-1130-4f76-a204-31b3ede61697",
   "metadata": {},
   "outputs": [],
   "source": [
    "miu1 = K_means(X, k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "66b7a95e-c441-42f0-b202-c82041065615",
   "metadata": {},
   "outputs": [],
   "source": [
    "Y1 = KE(X, miu1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "5692e4cc-9651-41b3-b86a-95c5eee939b9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAwD0lEQVR4nO2df3Ac53nfv8/dgZTx60CDNMmDDeIXndT1mBAIk5TGJtiJ0yqaRJp46iYtkypNZ1grNkGKTDNxOTXJuJzpTCL+dDyJpvY0tdg47dgeK6kmsT0TQnJNUgIIyLUtmyQAihJAiiJNgOCREnF3b/9YvIvdvd273dvfe89n5oa6w/54dwV897nv+zzPS0IIMAzDMPElFfYAGIZhGHewkDMMw8QcFnKGYZiYw0LOMAwTc1jIGYZhYk4mjJOuXr1adHV1hXFqhmGY2DI2NnZTCLHG+HkoQt7V1YXR0dEwTs0wDBNbiOgNs8/ZWmEYhok5LOQMwzAxh4WcYRgm5rCQMwzDxBwWcoZhmJjDQs4wTLI4fRro6gJSKeXf06fDHpHvhJJ+yDAM4wunTwO7dgH37inv33hDeQ8AO3eGNy6f4YicYZjkcODAsohL7t1TPk8wLOQMwySHq1edfZ4QWMgZhkkOnZ3OPk8ILOQMwySHI0eAxkb9Z42NyucJhoWcYZjksHMn8NxzwIYNAJHy73PPJXqiE+CsFYZhksbOnYkXbiMckTMMw8QcFnKGYZiYw0LOMAwTc1jIGYZhYg4LOcOEQR32A2H8g7NWGCZo6rQfCOMfHJEzTNDUaT8Qxj9YyBkmaOq0HwjjHyzkDBM0ddoPhPEPFnKGCZo67QfC+AcLOcMETb30A+HMnMDgrBWGCYOk9wPhzJxA4YicYRjv4cycQGEhZxjGezgzJ1BYyBmGqR0rH5wzcwKFhZxhkkLQk4vSB3/jDUCIZR/89GnOzAkYFnKGSQKVRNXJMZw8CCr54PWSmRMRSAgR+EkHBwfF6Oho4OdlmMTS1aWIt5ENG4ArV6rvb8wyAZQIupL4plLKQ8MIEVAq2Rk14xAiGhNCDBo/54icYZKA28nFWrJM2AePDCzkDJME3IpqLQ8C9sEjAws5wyQBt6Jq90Gg9dEPHACeeop98AjAQs4wWuJaVu52ctHOg8BsQvWv/krZplRSvHgW8VDgyU6GkdQy4ZckTp9WouyrV5VI/MgR/XW7nVBlXOPbZCcRfYiI/pGIfkpEPyGiPW6PybgkrlFl2NR7WfnOnYogW0XXXK0ZWbywVgoA9gshPgJgG4DPEdFHPDguUwte5BPXK0ELVRAPXLfn0O6fspALzlIJHyGEpy8A3wHwq5W22bx5s2B8YsMGIRQJ1782bAh7ZNEnyHv3/PNCNDbqz9PYqHwelXOY7W98eT1mpiIARoWZ7pp9WOsLQBeAqwBaTX62C8AogNHOzs5grjrKPP+8IhBEyr9e/TEQmf/BEXlz/CQThLhKgnhouD2H1f7ptPe/t4wtfBdyAM0AxgB8utq2dR+R+ykYHJG7w68HrJEgHrhuz8FBQeSwEnJP0g+JqAHANwGcFkJ8y4tjJho/J9W4SMMd1Sb8vCKIqkg756jkoXPlZmzwImuFAHwVwOtCiKPuh1QH+Dmpxs2K4kEQD9xq56g2Mc5BQXwwC9OdvAB8AoAA8CMAE0uvxyvtU/fWCtsf4RKUfRKFcVQ6h53fw6jcK0YIYW2tcEGQQ3bsUP49c8bFQeq98CRM+N4vw90LYwd3P4wSbH+Eh935iXooqmIPPDGwkNtkxw7lNTKivOT7mglqUo3Ri7JZiTmgn5+w8o7/4A/CF3cvHzD16oEn8SFt5rf4/YqjRz40pLykjSjfMxHHTlGL0Re28o6N6XhBF8P4kbbqlQceFy89yFoBH0AQBUF2X3EUcgkLeMywEuVKf8hW+dNhTFBrBTKdjsYkuVG0n346PuIY80QDKyFna4VJNpVSOq3mJ5x4xH42jDJaPMWiuzF4YSmY2U5/8RfxaTaW0MZfmbAHEDdcZaswwdPZ6bz16pEj5ZktROYZHn5ODJpNzJphZwzGbB3p+wPO5mfMxmR2X4BoiqPV70PMJ3g5ImeSzeOPKyKspdqEnllW0Wc/G/zEoB0htDsGr6qJnYhzFMUxqRO8Zn6L3684e+RMhLHj3RIpn3txfL89YC+bVnnVNyWqE8FOiMvErAngyU4m0ZhlI1iJV0wmtjzNsPBqks9qTE8/HVtxjBNWQs7WChMvrCbs4u7dmmFm8Tz1lHKtTicsvbIUrIrZvvIVrosIEzN19/vFETlTE5Ui1CilDJqN26tcbbcLRQQdNcfYxogiYGuFiT2V7IGoerdRtEeCIubFN1HESsi5aRYTHyo1efr6182bYT31FPDii9Yrw/uNlyvPx63JlZfXzgDgpllMEqjU5Cmq3q0XBShyXsAq6JILI5t55mH2FXF77UnsieIXZmG63y+2VhJGUD5oHL+qu7VD7PaKMbsfYd8vN9fO2TGmgD1yxheCFou4TZ5VSousNH55nXYF3Ewow/bU3fxuRHXOI2RYyBl/CFss4oBWlO0IkdMo3EzshIjG4sm1PnijnIUUIlZCzh454w6vmxAl0ReVvec3bCj3uc3K5O32WLFCziVYzSlU8tS9pta++1FpXBYTWMgZd3i5yky1xYDjjt2Hnhth0hb5mBUBAUoXxajfX7OxG3vmSKLY0yVgWMgZd3jZhMirxk5Rxe5Dz6kwpdPmLXmNmTzpdPm+Ub2/UWlcFhfM/Ba/X+yRJwyvJiCj4Ola4cU12p38c+KRO5nsi/L9tUvcJrs9BjzZyQRKLX9wUZ04rSX7wur67d6X558Xor29+qSmk06Otd7fOhfPKMFCzgRHrWln1fYLS1CcCqCXKZnV0hCdPORqfSDFLXc/wbCQx5DYrg/qthDEKpL1UlCcPBScWhJ+fLPwyhYxRvrt7ZWvParfkuoUFvIYElsh98OL9VJQnD4UnJ47ytfv9NqT4KtHBQ++UbKQxwgp4PJvJnaCXk10avmFrlQg4rcX75Xwt7dXH5sVXn0jcXrtHJF7g0f//1jIY0TshbzSL22tv9BelmzXEmU6efg8/7wQK1aUH7+hwZ237MUcgdNr5xYM3uDRA5GFPIbETsC1WP1Busmc8GoptyCiTKuMk7Aj2VqunZuiuccji4qFPIbEWsitcPMLbRQUK6ul2rGCEIyoestRFssk2zg+R+Rc2RlhzpxRXonCy5L+9vbajmXVu9zLXuVeXmc1nPSnCeLaa8Xrvj1RwssKaDPM1N3vF0fkdYyXOeYNDeVetPFYYXmuQUW+Ts4Tdf85yRG5EJy1wiQML6s+29utj+WnmNq5hiCE0674RdlSkcRhjCHjq5AD+BqAGwB+bGd7FnLGMbV4zn5FeFESHLv3JS7RbtS/NYSMlZB7svgyEW0HcBfA/xBCfLTa9rz4MuOYWhby9Wux4igtKmx3LHFbuJkxxdfFl4UQLwH4hRfHYuoQO5N1tUwW+TXhWGlSzs3CGLXsa/e+BDn5GhRBLkIS9QVPzML0Wl4AusDWCuOUWifr2tuVVzWP2g8LpJJfX+v53IzVrl8fFTvIC4K8ngjdO/g92VlNyAHsAjAKYLSzszOQi2ZiQK0FKn5nalTaz+r8ToqAjMcPooAoSf5zkJ5/hOYXQhdy7YsjckallpLxdNrfPyw7DwozUbR7LWbHt3qFXUAUVYIsuIpQcZeVkHNBEGPJjh3Ky1eceLdyTc9i0XwfrwpH7Cw5Z7aosN1rcbK4cpw9bD+xc6+98rVjML/giZAT0V8DOAvgl4joLSL6914cl6kDnExiVhNAr/6waq0wtHstTh44jz9uf1s7RH3Szi7V7rWXC3n7XZXpBWZhut8vtlaiTeDdF+16t5Va2Xo5+WTliabT1ufQruQjrR+ra6nUJ8Zvjzwik3aeUOn3xmtfOyLzC+DKTsYuToU8sOZetQhsLVTysO0ullxJIN145G4EJahJuyBFz+pcEfK1vYSFnFGxK7xm29n9zBeCTjmzO6laa+aN1fGtjuH2+quJmxcCHJW0wAhlmngJCzmj4pWQh7IARpDRnt2ortbor5JVZCaAbsSp2oPJKwH2U0CdpGwmzUZagoWccSW8ZvtmszFfyagadkXJznZmD6BKXrnZQ6rWB4Ydq8grAfbL0qjFjoqIr+0lLOSM50JujMwjh/EP+emnnf1h243qqm1n9fOnn/Z3Eehq+2nnFrwSYL8i8rAmiCMGCzmj4kZ4Q/XInWAngrPbB92O+NeaQWG1n9nntdoFdkS6ktg7iWj9sjSq2VBenivCsJAzKl4Ludd4cg67EVwQ0ZuXCx7XYhfYtX68ePDJY3ltadTSjz6BsJAzscETIbcbwQWRjubUbnC6fTXhdGIRyeP43QbBKQmdvHRKIoQ8kl/hGc/wNAsmShG5UxFyEsHXItJ2Itco5mEncPLSKVZCzr1WmGRy5IiyaEIlgiqzdrrgsZPeHnb6wsgxGHvDVCKK/UWcXkM9Yabufr+cRuSh5CsnFK/vnR//Lzw7ZrVIPOiIzsnEqd0IPsh0vzq0MqIGOCKPL4F0IUwiGzZYfx50RGeniZNsaPW7vwu8731Ae3v1CL5a5Fxrkyyn3yKYcDFTd79f7JE7w4vr9vpbTSy+JUUpqqyUdeFmrNUyXKJy/YwngCPy+CEj8ZER5cWRuUOiFFVata69dUuJku163UYqXWOtx0wQivZZv08KmbAH4IQzZ8IeQXyR904+CNzeS6+P5xs7d+qFW1oNV68q9sORI8EIe2en+Wr3gCKstfZAr4Qfx4wR09OHUCjMoa/vGIgIQghcvvwMMpk2dHcfCnt4nhIrIa83YiOWMeFL/+Q0/uPFXXiotBSlSp8a8F/MjxwBfud3zH8mHypmQl8tS0R67/dMrqnWYyYAIQQKhTnMzJwAAPT1HcPly89gZuYEOjr2QAgBqpbVFCMojK8ag4ODYnR0NPDzxhUWcm+4/lAX1r1nImxy8tNvVq9WrBSz8x85ohdkQEmPrGYFdXWZi7WbYyYEGYFLMQeAjo49aoQeR4hoTAgxaPycPfIYcOYMi7gb5NzCB94L2Wo4ccJ6ybBa/fxK9kmU5ghCgIjQ13dM91mcRbwSLORM3XAVIRe5VBPWWgperMaeSimvAweUB0UdFtHIiFzL5cvPJHLCk4WcSTzyG83Xf/kI3k2FvIiu19WJZgsDA0CxaJ2vXgtO89FDXuRZa6t0dOzB0FAJHR17MDNzIplibpaT6PeLm2YxYTA0JMSf/HIC+3X43ezKizVJQ8hfn5o6KC5e3CNKpZIQQohSqSQuXtwjpqYOBjoOL4FFHjlPdjJ1Qd1MGKdSinQaIVK+BdRCpQlVwySxEALU3W17e78RhuwU4/u4YTXZyemHDJMk/Eg5tJmPruZtX70KU6kMIX/dKNpxFvFKsEfOJJq6q44188zdzgPY6IQoNHnbhfXNzo7DuIaFnKk7JiYSLOZ+pBzaeDjIVL+Ojj249O8WUFyJitsz3sIeOVMXaD3yuvHLvUT2g6nS2kAIgZGRFD7wfaDnvwErbxAoyFYICYc9cqbukZH4yIjyngXdAcaeNSYITd72jU8pr46O4cQW4UQJtlaYuuDMGaC/P+xRJBcp4nWTtx0xOCJnYoEX0TM3IfMPIkIm06brZSLL4zOZtshF5ElLS2QhZxjGE7q7D+kEUYp51AQyie1tWciZSCOjZy99bT8j8XqP9qOet61NkwSS096WhZxhaiAugp00C8EtWstnZuaEKuhxb2/LQs5Emrj42n58c3BLEi0EL5Biru1THmcRBzzKWiGix4jo50R0mYj+2ItjMkwUiUulqNZCkFkj0kIoFObqOotEmyYpiXtmjeuInIjSAP4cwK8CeAvAq0T0ghDip26PzTCSqEbikqh9c0iqheAWY5qk1iMH4huZe2GtbAFwWQgxBQBE9A0ATwJgIWcSR9QEuxJJtBDcErc0Sbt4IeQdAN7UvH8LwFbjRkS0C8AuAOjk5jlMQomSsFtZCPUu5nFJk3RCYJWdQojnhBCDQojBNWvWBHVahvGFqK+jypWWlYl6mqRTvIjIZwB8SPP+g0ufMQwTEkm1EBhzvBDyVwFsJKJuKAL+2wD+jQfHZRjGBWFbCJzDHhyurRUhRAHA5wH8A4DXAfwvIcRP3B6XYRj3+G0hGC0a+X56+pDOwpFWz/T0IU/Pzyh4UhAkhHgRwIteHIthmHhgXXCURaEwn7gy+CjDlZ1MXRCHdMEoYxTfUqlUsWdJb+9RAJzDHhTcj3yJqFboMcnEypKIImY2yeTkPqTTWTUTZmQkpSuySaVS6uSqhEXcP1jIHRInwY/TWP0iiiX1cfKPK5X6F4vzauQtMdosWjjt0T/q3lqJYrMjJrnErY1qpVL/3t6jmJzcp9v+8uVn1M+TVgYfZepeyO0SJ8GP01j9Jmol9bX2QAkzlc+s1L+aWEvbhXPYg6HuhTxqf+hMsITx/91JDxQhBK5cOaxmh0iCbEdrZpNMTu5DJmMt1nZz2DnX3BvqXsjtEifBj9NYgyJK98BuD5Tp6UNYXLwNImBm5iQAASGAhYXzWFg4H4gVU6lboLRXrMS6Wg4790v3DhbyJaL0h874T1j2k902qtJLn509iVxuGB0dw0tirtDRMRyI11yt1D+VSpVtb4eg5grqJeJnIXdInAQ/TmOtF7TCKKPZZWHMlkW3AHQWjKSv73hgguRHqX8Q/dLrKeLn9EOmLpHdC4eGlFeQ3Qy7uw8hnc5icnKfKpC9vUdRKMzrUhC1Ymck6FQ+K5vETT682fVpRdcNTlZIilNOvxUckTNMwAghUCzqS9i1GSBS3IUQuHRpr27fXG54yTMPP5XPbcRrNVcgM2LcRM52I/6kRO0ckTN1TRh9xaXIWFVFagVldvYkWlq2oqNjGLncMGZnT0IIxSPXWjFAsJGk2zVBtds3NfUDAJqa+jEzcwJjY5uXjnNbVzTl9PqqRfzaa7h0aW+s1zXliJyJJMbqy6T5/dVSEI2TjNr9Mpm2JRGa10XvQUaS0g4SQugi3lzO3iQsEeHu3Qk0NfVjYGAUU1P71WPk8xNoaFinE/FLl/ZiYeE83v/+x2xfn1XEL5t69fUdQ1/fMQghMDt7ErOzymRyHHvCcETOMCFgp4S9u/uQKijy1dd3DF1dB1VrppZo2Aumpw9hfPxRGLVuYeE8rlw5XHV/IQSam/uRz09gamp/Wan/4uJ1zM6ewuXLe3Hp0l7Mzp7EwsJ5XZRe7fhWKyTdvPmCeu8AlF1D3EQc4IiciRjGtEDj50mIzJ2s5G41yWjH//Ur9U4IgcXF22o+u5aFhfNobd1a9VxWHrakqUkR+fKUS3vZOpXSJtPpLNrahkzPC8RzXVOOyBlbRKHZVFKwEpmOjj22S9gr+b+A88ZclTI3zCLgvr5jaGzcZHEse369tGe0bN9eREfHHuTzE2Xb9/Qc1WW0VDuH9huNPF9f3zH09Bwuu3e53HCs1zXliJyJFMaqVOPnScFtbnal6lAAjoptpqcPqp4xEaFUKqlZI/JYZlkdq1c/iatXX9ONIZcbRibTpu7f1XXQ8luBbIerZXJyH3p6nsWNG9/A4uLbup+dPbsOa9f+BwB3kU63oli8YzonYLyvZu+N904OMa49YVjImYpwAy7/qHUZNjvWjN1im6mpg7h16wU1Au7tPYqxsc3I5yfQ0TEMIaBOAuofCMMolcwjVlmR2tKyFYXCbfT2Kv3J5bjT6Sy6uw+VlfrLFMy5uTOqiBM9BCHeXTruTczM/FcARWQyq1Eo3Cx7MFVLJzS/d3uXLBxS712cRBxgIWciCj8orKlcNr+cklitMZfMZ8/nJ9TUP7l9U1M/enuXJ1r1DwRFxK9dO6W+v317BPfuvaaK/vr1u5FKEWZmTmJu7iVs3jymCnVTUz+IYHkNd+9OoKVlC+7fn0KhcNNw9UUAiqg3NfWjp+dZXYS/uHjb4sGzLPjl5z0OgGIXhWuhMHygwcFBMTo6Gvh5k0aQ0TFH4tHDaJFMTSnZLFIQlyNNBbOIXBuhatm+vaj2URFCYGRkeTptaKiEK1cO4xe/+Hu0tm5FX99xCCHw0ktp3TZCCDW6l8hJTG2LAjPrRQiBUqmEl1+2jjU7O/8zisU7huh7L+7c0U/CWl13HHuwENGYEGLQ+DlPdjJMTDGLrmVxi1bEl5tulU/imU04AlDbB1h58V1dBzEwcHYpmoXpAhNEhM2bx3Sfa0V8cnIfrlw5ZHl9xvMauXXrb01SME+itXWrbrtarJK4le2ztRJDwvCtORKPNlZNtnK5YWzceFxuVWYflEoljI1t1h1L2iyybe7s7ElLLx6ApV+vRLnlY+3peVZns2zYcFDjoe9FJtOGa9e+igcPZpBOt6NYvFV2jExmtaklJH19LcZ0wmo+ehzL9jkiZ2IFp0FaY5aSuHHjcV0xkVaIZNaIFERt6l9TUz8ymTY0NKyqmCZp5dfncsNYWDiPmZmTagm+5OWXM5iZOYGGhnXI5ycwOfmMKuIzMydx69bfQ4gCAKBYvIV0ul23fzrdrnrk7e1P6H6mffCYpRMqFbG31c9KpdLSeZViqmKx6Kr1QFiwRx5jduwAJiaA/v76iZjZq7fGzO+uVm4uo8/e3qNlmSU9PYfV41bzk822mZ4+pGbEdHTsQU/PszrPO5NZj0LhWtmYGhpy2LbtKs6d+yAWF6+rnzc2bkI2ux0NDW0oFOaRTmdRKt3RXW9z8xa0tm5DX5+SKaNNpZTRtn6xDoWWlq14+OEfqisfFQrzju5jUFh55GytMLGA0yAr46RaVIudfHY7aZJm2/T0HAYR0NY2ZLpQM7Boei0rVnwARIRHHpnRTaAODl4AEeHKFeW4xeI8ZmdPorl5C4R4gMXFG7h79xW0tm5VrZBicR6ZTFa9TpkauX79bt05m5u3YHJS8dilh29nKb6owEIeU6SQzc8r4hZ1YYv6+OJO5ZTEyml11YRam0li9m8lursPq1Gx9MSldaNks6QhUwqVsa7GwMAohBC4cEEfeI6NbcbAwOiS9aHkqedyuzE//zLy+QlkMqsBAO+887+xuHhdlyEjx6ptkqVlOZVyj+lDJ+pl+yzkjG94Kd68Dml1/FjJR1ov6XQWxeK8KnLGyspK9ksqlUI6rSzU3N39Z5ie/kNNtFvUna9QuInx8Ufw7rvTqg8+MDCKs2cVP/3ChUEMDCi27MzMCTXNcPnBANWOkSJe/g3D+nq1E7FOvtmEDQt5TImLsE1MKGNkSyQYaq0WNUPbr1sK5dzciBr9yspKMx/a6LsXi/N4551vY25uBAMDo6bNqiTvvvsGCoWbanQ+NbVffZ9OtyKdTpcVO23ePKazYSRmwnvnzvmy7SRTU/vVh47TbzZhwkLOeI7Rz56Y8O7Y/AAIFmNKo4x6ZbSsLa2Xoi7Ff25uRFfRSfQQHjy4irNn1xnOQli//nOYn38Z9+4pvVuamjYhn39NnRzt6BhGT89RpNNpFItFTE3prQ9jCqVEa4nIeYSFhfO6ylNAW4l6QlesBHjzzcZvWMhjThyErb+//rJr4o42l9oY/Ury+Qk1CtZGsL29R9XIXf5cRvREDxnK7lMASrh27ctobt6CxsZNKBRuIp/XN+MqlQRSqRSmpg5idvYrSxH6Wmzb9hbOnVu/9IBZCeA9nQ9vtES08whXrhxW884bGlahq+sgZK69rGqVRFnEARZyxgfMbB/O/Y4PWkvFqqjHiDZiTaVSZVZHPj+xtN5oC2ZmjqifP/rofUxMbMG9e6/h7t1X1M+1njegTEYKIbCw8AP1QVAovI0LFwaX3qcAvIeGhrV4+OFXceHCINLp1rLWwIpYK8Lc3X0IpVJJ1yYg6pG3FSzkTCBwJB4cbvuIVMrukAIrPXKJmYVhNq65ub/TffbDH64EoOSIS1sFgNp9sVQSuHlTyUK5fv3LZcdc3qcEAFiz5rcwObkP9+69hpaWrapPD+i/ZcjxGBd5jqOIAy4rO4noM0T0EyIqEVFZkjoTP7ysnAxjYeN6x+6CEla9ROS/RKQp7Vfo6BjGwMAompr6USjcREPDOnR2flFXPVksFnHhwjbMzp4sW3ji2rVTOrHWYvy8qalf7bC4evVnbF//7OxJNZWwtXWrrjOiHxWbUenJ4jYi/zGATwP4Sw/GwjCMC7RiBVi3cbXqJXL37gSam/vVib7Ll/fqjn/79siSDbECANDe/i9VqwIAFhbGMT7+KN599woARZwbGzfhwYPrKBT0i0QYI3BA6QsDKGKcz08gn59QJyGd0tExrLbhBayXlnNTsRmlniyuhFwI8ToQ368jzDJcORl/7IhVJbGXk4NzcyPIZj+J2Vklsm1s3ITFxeu4d+811fdubt6CfP4VXLjwCLLZbejufhbj44PqJCVRE4TIq2KdTq9GsbhsxVhF5n19x3R2jjazpKnp4wAWTZeBMyJFtaFhlS7PvVqPdrvYfWgGBXvkTF1QLw+mamKlbVtrFPuenmdx4cKgGg0D2v7hw7reJESkFuPcvfsK5uZe0mWaCJHXjatYvImOjmGk022Ynf3Lsgg9k1mLTCZb5q3Pz59DKtWEUimP1tZtyGTacP/+VZRKvwCgpA3euPE/yzokyodQLjes9mm3qtjU+ujK2KuLsB8RvhuqeuRE9H0i+rHJ60knJyKiXUQ0SkSj77zzTu0jZnxB+tlDQ8or6f52UrsoWvUPl97t9PQhTE7uK+tBnk5nkU6nLfqHl7eGvX9/umy7apRKJSwu3lZFvKFhrdrzpFB4G7du/a3aO0V2TLx79xWUSspD4datb6FYnMe6dTvVY167dkoV8XR6te58DQ3r0Nt7VO3TPja2Wf3mASy36x0b24xSSZksrbZItRazbpNhZb1UFXIhxKeEEB81eX3HyYmEEM8JIQaFEINr1qypfcQM4wAp2CMjyqutTXklEWPjLGMbV22xzujogG7fmze/g2KxaNLYCmr2imx1m8sNo1C4YXtcctLz2rUvqxORALC4+DauXTulLnyRyWSRyw0jm91m+mB48GAGs7MnMTt7SvXTJevXfx6PPnpd91kut0utApXteYHl0n05cau00t3neAK02kMzSDxpY0tEZwD8oRDCVm9abmPLAMHYHUbvP72U2lxcavExNOT/GIKk2gRcqVTC2bM53Qr1cuIxk/kACoUbZTncmcxarFixDvfuvaZWPb766sO4f/9HluNobNyEtrbtOp/dyhffvr2o5nLLJd7MWtha7Q8olaBEK3W56Ma5AeNydUSka+hltp8VlbpN+mmv+LLUGxH9JhG9BeARAP+HiP7BzfGYZODWtvAjBTKbVV7F4rKIJ5Hu7kNlnrh2QQkiwkMPden2WbVKeZqVSnlkMmvVqshPfrKwlGqoiH5Ly1bMzJzASy+lK4o4IFvXHkNTUz9SqWY8eHDdctuxsc269McLFwZ1Ig5YT47KMebzSkHR+vW7Tb+JWEXOqVSqJnvEqtuksQApKNxmrXwbwLc9GgsTMkFNCEYhQyabDf6cXmNV+FOtcVZz8xbd4sRyFR9FwP8pikWl0+DU1H4MDIyqE6By1Z9qNDSsRSbThqmp/ZriobfR1NSP1tZPqvZKY+MmCCGQz0+oLWr1+9wsO+6aNf9KjfIBqL1YGhrWYeXKTnz4wycMja6yFbsZumlZ60e3yVrhrJWEEoY4uhVop/s7Of7cnH6fuFNLDrPc5tq1U1i/frfOr5aTmn19x1EqlTA1tV+XjaGU15cfM5fbjd7eY5ia2oeZmZNIpRqxuPg2rl79EwBYishbcP/+67psmIaGtUinV2LTpv+L8fGPI5+fUEVZPlRyuWFdKmJDwzqdP6792bZtbyGVSpmK6vT0IdNuhul0ZZG3G5lXIqg0RF6zMyG4sSOME4I7digTgn6JXhQyZOKclVNrlaK0A3K5YdMim95eReCmpvYjk8ka9lUi94aGdcjldmPFihwApUhofPxRELUik1mNUumebr98fgKl0kJZdL1mzW9hYeEVnD//Ibz//b+h+1l7+2+YPjiUtrnDZVaIHLPxHkkB7e4+ZNrNsKfnsKf2iN2qWj/giDxhhGlbuO2Rbnd/N9cYV/HW4iaHuavroLrIsXFSc2xss2p9GBdMvnPnnDrhODt7Cg0NawFA9cpLpfdQKNxEKtVoKuZGq2RuTvmft7h4HW+++SXd9oXCvPrgWI6UlTGXSgKTk8+oCyxr2+gCULsaGr+tWPVU8coeCbtAiIU85hhFLZu13NQSrYDK3uFBLCGXBFENi1qrFIkI6XSbrtBH9jSR1kdj4yb139Wrn1TzsHO53chmt+PatVO6rBdgeTLSKOISpX/5JrVoyGzycvv2oirKLS1b1ehbKWI6Btlidm7ujNoLPZVKqW1z796dWDpXZUGV96HSPXJK2AVCLOQJo79f/74WsXS7EIRbga62f1xWR/ITqxxmO6JhXBRZCKHzy9vbnwCgiO3t2w+hre2fo6mpf6ln9yHdtlqamz8OIcxL6Jua+tHS8omyPuNatIVKmUwburoO6nrDLPeAuY2ZmZMYH38UAwNncfnyM2puOADLytW+vmOYnj6EYnHel/4oXrYAcAoLecyxErVa/G1jZM4LQUSTSjnMgD3x6O4+jKmpgxgffxTF4ru6n2mtjtbWrSgW55HPTyCb3Y5Ll/ZUGNeibmFlrXWjnejUItfklJOrAHQl82aWhZwCWFg4r+aG53LL3rlio+i/msq2vLduvaCOw2v7w83D1S0s5AmlFgE22jRyvc2oinlUx+U3VjnMgL11JaW9UCzOqemEjY2bQARdxJzL7UZf33H1vTZLJJ1uL+tvIsW7vf0JNdqfnHxG16MFUDoTzs+fw927ryylOe5TJ1qNq/NYWRbGjBZ5ydoGYFouXdq7dH3LKwd5aX948XB1Awt5QvBD1Iw2DRMdtJN02s5+2vdm6JdwO447d85jYeG8RcHN8gTgxo3HdcLZ2NiHhYVbaGzchFRqpdpEi2iFzqIwJtAo3rYyzosX9+Du3VeQyaxSC3PMxm1mWRg3m5k5qemSuJy+KCdN5dhle1vt6kVeiKzbh6tbWMgZFfae40Wl3uJmnq+5TWHdokOKX1/fMYyPP1p2rKamTWhvfwJdXYfU/uWZzCrdOGR2SV/fMVy6tBezsycxObkPfX3H8OEPn1CvQ/uvPL52UYhLl/bqzq9ktCh57/K4EvmNYDmSX/5Zb++xmguAqhFmgRALOcPEFKcpb1Y2hRkyql1YOI9Ll/aqFkxLyxaUSg/UnibZ7HZMTj6DTGYV+vqO60TMGKFu3Hhc/bySuBknN6VQNzdvwebN53DhwiNYWDgPIaTg6/eXLWsBlHnWY2Ob1YlRP+yPalW1fsFCzpTBkXg8qCXlzcymMKOl5RNoa9uOTGYVALl6j9CVxzc3bwERqfneRmqJULUPp7m5EQwMjKoPESEeYGrqIAYGzuLSpb3IZNqW0hWXo/NlYRYQArpvBFr/XD4kgrQ//MST7odO4e6HDOMdVp39rLbVRqEA8L73fUzXBKu5eQva239NTf+T+wHQnUciC3O0k5RuKJVKauQs0WbCrF69bOdUspYAWPwsi+7uw7p7EhcRt+p+yBE5w8QYJylvZpkVFy/uKcsLz2a36URce1wzZHVltVxs4+Ss8V9JKpXCww+/ih/8oEH9TJvO2NY2pH5eKerXFv/I/zabEI6LiFeChZxhYorTlDejbw1A7bnS0rJVLa5R9qeySHZm5gTWr9+NO3deLivtl76zVXQrI+dMJotCYXnZtXQ6i2JxXvcQmJo6iFu3Xig7htbbNl6X9p4Yi4iUkv3bEAJLRU0HQ1sk2S9YyBlHcEZLdKgl5c0YwWYyq1R/2Wp/bbMtbS62tthH6zsb0frecr+5uRHdceRDQFu0Y9bKtrv7z6qmVvb2HtX57KnUCiwsKJOzudyw2rcljEWS/YKFnGFiTC0Titqf2d1fbnflymHVE9fmYre3P2HpkRsnZQHoHgLaBxERob39Cbz33ltlIg4A4+Mfx+bNY+q2EmMGj+y/YqwmXc4nD2eRZL9gIWdsYaz65Mg8OrhNeXOyv1wuzpiLXSzOV4xuK2XLGAW1p+cwhCjhzTf/i/pZLrcbc3MvIZ+fwOXLz4BI+Tah7WZoN7XS7Jxxh/uRMwxTFdlrW7vGZVNTPzo7v6hbVs0qC85sUlZi3E9Zt3NBt83s7Cm0tQ2pqxTNzJws672uFfNq2Fkk2fjzMDL87MIROWMLrvqsX4y2RSaT1WWQKMU3QuerGysztTncxsZa2slZALoJ3N7eo2qFqLZ608waMXtYNDX1I5v9pJr/vly2X7kIqJYVmMKEhZzxHRb/aGO19qfEzOMGlsX0ypXDEEJZtELurxU9ZbI0u7TyThZtbUPo6XkWU1P7kU63oq1tSH0ICCFMK0IBfdMuKxGXD4B0OqtOmqZSK5HL7QZAataK7G1udt1hLxJRCyzkjCNYjJOF3cjTzOOW4l4ozGF29qS6jVH0lPS/ebVoSNozCwvjaGl5WP1ce+7yFEP9uI258mYZPF1dB9X8dvmQkdsbj2+8D3KyNIxFImqBhZzxDZ4gjTZOIs9KhUeV2gQAyyv2CCGwceNx1WNvbNxkeW7tGJXPy0vt5X5yjMYMnGodFSvdh8nJfWUZL1EVcYCFnGHqFru9WuwUHlVaGaev7xju3Dmn87mlV51OP1Tx3E5z5WvJ4LG6D8Y1TYNaJKIWWMgZ3+AJ0uhjZZlUqgg1iilQXr5vFL3W1q1q8yvlmMutaLWfV8phd5Ir7xSz++Bnl0SvYSFnmDrGbq8WKzGV21eK1pXz6M8rRdz4uVXU6zZXvhpWGS9x6ZLIeeSM75w5w9F4FDFaJkNDpYo54WZiahWtKxkqbQCgpg8qrXCXuX17RG0zW+3cfmJ1H/L5CUxO7tOtwBTF1EOAI3KGqVtq6dViRjXrQ9unRUuhcFNdNDnMqNfufYhiJC7hfuQMU+dUyyP34vjaRlWK/bIX2gUhzAqJgsbv++AFxP3IGYYxw2//WYl4Vxki3uMwK8oJUzj9vg9+wkLOMIzvhLkwcT3Ak50MwwRCnCPeqMNCzjAME3NYyBmGYWKOKyEnoj8lop8R0Y+I6NtE1ObRuBiG8ZA49dZmnOM2Iv8egI8KIT4G4CKAL7gfEsMwXiIXhZDiLQtgpqcPhTswxjNcCbkQ4rtCiMLS23MAPuh+SAzDeIW2s58Uc1nFaFxhh4kvXqYf/j6Av7H6IRHtArALADo7Oz08LcMwVtjtcMjEm6qVnUT0fQDrTH50QAjxnaVtDgAYBPBpYeMRz5WdDBMsQgiMjCx/AR8aKrGIx5CaKzuFEJ+qcuDfA/DrAH7FjogzDBMsdjscMvHFbdbKYwD+CMATQoh73gyJYRivcNrhkIknbj3yLwNYCeB7S0/2c0KIz7oeFcMwnuBVh0Mm2rgSciFEn1cDYRjGH7jPSfLhyk6GqQO4z0myYSFnGIaJOSzkDMMwMYeFnGEYJuawkDMMw8ScUNbsJKJ3ALyx9HY1gJuBDyI68PXX7/XX87UDfP21XP8GIcQa44ehCLluAESjZiWn9QJff/1efz1fO8DX7+X1s7XCMAwTc1jIGYZhYk4UhPy5sAcQMnz99Us9XzvA1+/Z9YfukTMMwzDuiEJEzjAMw7iAhZxhGCbmRELIiehLRPQjIpogou8SUS7sMQUJEf0pEf1s6R58m4jawh5TUBDRZ4joJ0RUIqK6SUUjoseI6OdEdJmI/jjs8QQJEX2NiG4Q0Y/DHkvQENGHiOgfieinS7/3e7w4biSEHMCfCiE+JoToB/B3AL4Y8niC5nsAPiqE+BiAiwC+EPJ4guTHAD4N4KWwBxIURJQG8OcAfg3ARwD8ayL6SLijCpT/DuCxsAcREgUA+4UQHwGwDcDnvPh/HwkhF0Lc0bxtAlBXM7BCiO8KIQpLb88B+GCY4wkSIcTrQoifhz2OgNkC4LIQYkoI8QDANwA8GfKYAkMI8RKAX4Q9jjAQQlwTQlxY+u8FAK8D6HB7XLcrBHkGER0B8G8BzAP4ZyEPJ0x+H8DfhD0Ixlc6ALypef8WgK0hjYUJCSLqAvAwgPNujxWYkBPR9wGsM/nRASHEd4QQBwAcIKIvAPg8gINBjS0Iql3/0jYHoHz1Oh3k2PzGzrUzTD1BRM0Avglgr8GRqInAhFwI8Smbm54G8CISJuTVrp+Ifg/ArwP4FZGw5H4H/+/rhRkAH9K8/+DSZ0wdQEQNUET8tBDiW14cMxIeORFt1Lx9EsDPwhpLGBDRYwD+CMATQoh7YY+H8Z1XAWwkom4iWgHgtwG8EPKYmAAgZY29rwJ4XQhx1LPjRiH4I6JvAvglACUo7W0/K4SomwiFiC4DWAng1tJH54QQnw1xSIFBRL8J4BSANQDmAEwIIf5FqIMKACJ6HMBxAGkAXxNCHAl3RMFBRH8NYAeUNq5vAzgohPhqqIMKCCL6BICXAfw/KHoHAP9JCPGiq+NGQcgZhmGY2omEtcIwDMPUDgs5wzBMzGEhZxiGiTks5AzDMDGHhZxhGCbmsJAzDMPEHBZyhmGYmPP/AXnPKYxlZHONAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "draw(X, Y1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "b39dbb53-1562-4a44-8b69-5518e2ca70f5",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████████████████████████████████████████████████████████████████████████████| 50/50 [00:05<00:00,  9.96it/s]\n"
     ]
    }
   ],
   "source": [
    "alpha, miu, sigma = GMM(X, k, 50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "f9d83627-21c4-4f26-945b-1a5f5a766b84",
   "metadata": {},
   "outputs": [],
   "source": [
    "Y2 = np.argmax(GMME(X, k, alpha, miu, sigma), axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "81945883-584a-46b6-81bc-b7ffa8baed17",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.2, 0.5]\n",
      "[0.38992172 0.49839234 0.11168593]\n"
     ]
    }
   ],
   "source": [
    "print(alpha_real)\n",
    "print(alpha)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "6428396a-4ef5-4f42-addf-0f64823e8734",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[array([-1, -1]), array([ 1, -1]), array([0, 1])]\n",
      "[[[ 0.05622333]\n",
      "  [ 1.07949114]]\n",
      "\n",
      " [[ 0.9961998 ]\n",
      "  [-1.01344529]]\n",
      "\n",
      " [[-1.28859858]\n",
      "  [-0.57686972]]]\n"
     ]
    }
   ],
   "source": [
    "print(miu_real)\n",
    "print(miu)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "84442b32-fd40-4e44-b966-bf3395392c5e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[array([[0.15401791, 0.        ],\n",
      "       [0.        , 0.15401791]]), array([[0.12400811, 0.        ],\n",
      "       [0.        , 0.12400811]]), array([[0.63992835, 0.        ],\n",
      "       [0.        , 0.63992835]])]\n",
      "[[[ 0.61872913  0.03126732]\n",
      "  [ 0.03126732  0.43276116]]\n",
      "\n",
      " [[ 0.11296431  0.0063979 ]\n",
      "  [ 0.0063979   0.12611463]]\n",
      "\n",
      " [[ 0.34984941 -0.18321799]\n",
      "  [-0.18321799  0.65368268]]]\n"
     ]
    }
   ],
   "source": [
    "print(sigma_real)\n",
    "print(sigma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "ea78a386-eb05-4ed9-8ffb-5a74b1a8bbb6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAvsElEQVR4nO2df2wc53nnv8/uUnb426FUS8uE4i/lAF8Q0RQjyUYi+tD0zglaGw0u1x9CL0UP8MVIRMlyUTQnXCRdz8ABqfXLaYAaSHC9Stf0DomRtGe0SYBadi6WbFKic0mdWCIpyV7KP6SIFLWUI+7ue38M3+HM7MzuzM47P/f5AAtpl7sz7wyX33nm+z7P85IQAgzDMExyyUQ9AIZhGMYfLOQMwzAJh4WcYRgm4bCQMwzDJBwWcoZhmISTi2Kn69evF/39/VHsmmEYJrFMTU1dE0JssL4eiZD39/djcnIyil0zDMMkFiK6bPc6WysMwzAJh4WcYRgm4bCQMwzDJBwWcoZhmITDQs4wDJNwWMgZhkkXp04B/f1AJqP9e+pU1CMKnEjSDxmGYQLh1CngsceA5WXt+eXL2nMA2L07unEFDEfkDMOkhwMH1kRcsrysvZ5iWMgZhkkPV654ez0lsJAzDJMe+vq8vZ4SWMgZhkkPTz0FtLaaX2tt1V5PMSzkDMOkh927gWefBTZvBoi0f599NtUTnQBnrTAMkzZ27069cFvhiJxhGCbhsJAzDMMkHBZyhmGYhMNCzjAMk3BYyBkmCpqwHwgTHJy1wjBh06T9QJjg4IicYcKmSfuBMMHBQs4wYdOk/UCY4GAhZ5iwadJ+IExwsJAzTNg0aT8QJjhYyBkmbJqlHwhn5oQGZ60wTBSkvR8IZ+aECkfkDMOohzNzQoWFnGEY9XBmTqiwkDMM0zhOPjhn5oQKCznDpIWwJxelD375MiDEmg9+6hRn5oQMCznDpIFaouplG14uBLV88GbJzIkJJIQIfadjY2NicnIy9P0yTGrp79fE28rmzcClS/U/b80yAbQIupb4ZjLaRcMKEVCpuBk14xEimhJCjFlf54icYdKA38nFRrJM2AePDSzkDJMG/IpqIxcC9sFjAws5w6QBv6Lq9kJg9NEPHAA+/3n2wWMACznDGElqWbnfyUU3FwK7CdW/+ivtPZWK5sWziEcCT3YyjKSRCb80ceqUFmVfuaJF4k89ZT5uvxOqjG8Cm+wkog8T0T8R0T8T0c+IaK/fbTI+SWpUGTXNXla+e7cmyE7RNVdrxhYV1koJwJNCiPsA7ATwRSK6T8F2mUZQkU/crIQtVGFccP3uw/j5jINccJZK9AghlD4AfBfAb9R6z7Zt2wQTEJs3C6FJuPmxeXPUI4s/YZ67kyeFaG0176e1VXs9Lvuw+7z1oXrMTE0ATAo73bV7sdEHgH4AVwB02vzsMQCTACb7+vrCOeo4c/KkJhBE2r+q/hiI7P/giNRsP82EIa6SMC4afvfh9PlsVv33lnFF4EIOoB3AFIDP1ntv00fkQQoGR+T+COoCayWMC67ffXBQEDuchFxJ+iERtQD4NoBTQojvqNhmqglyUo2LNPxRb8JPFWFURbrZRy0PnSs3E4OKrBUC8A0ArwshjvgfUhMQ5KQaNytKBmFccOvto97EOAcFycEuTPfyAPAJAALATwBMrz4+U+szTW+tsP0RLWHZJ3EYR619uPkexuVcMUKIEDxyL4/ECrmqL3WYk2qMGT73a7AHnjichJxL9N2iMj+b7Y/ocDs/0QxFVeyBpwc7dQ/6kciInO2Q5GK8k3LKhzZGoU5R++OPR28zqLQ6mvXuJMF2Edha8QnfhiYTN0Ut1guy00Xb+h0IW/SCEF6VdmESxDHhFy8Wcr9wRJ5MnH5vtf6Qa0XuYf/+jQKZzcbjO2gV7ccfT444JvzvmIXcLwm/kjct9ewUu+jRjfiHcUfm9m7C7RhURM12Y3I6x3EUx4TfWbOQqyApt4/MGo1EYHERK7cXFDdjUBWIxOUi1ygpjcg5a8ULYVX9Mer4zGe0zCAj9Ypa7LKKvvCF8Itj3BSIuR2DqmpiL0Vrccx+SWuRk526B/1IbETOxBs33i2R9rqK7Qd9R6ayaZUqSyGuE8FeSPCdNdhaYVJNXOwQlaicl1FlKcQ5NbMJcBJytlaYZOFUqGNnHQhhv42krGhjZ/F8/vPasXotVFJlKTgVs33962w7Romdugf94IicaYhaEWqcUgbtxh2H1g5RWAoJtjHiCNhaYRJPLXsgrt5tHO2RsOCUXeU4CTkJp9vPABkbGxOTk5Oh75dJOJmMvV1CBPz1X2u9b4z2SmurZkU8/7zzyvBBo3Ll+VrHX6k0MrpgUXnsDACAiKaEEGPW19kjZ5JDrSZPcfVuVfSel/MCTkGXXBjZzjOPsvmX32NvhsZlqrAL04N+sLWSMsLyQZN4q+7XDnFb3Wl3PqI+X36OnbNjbAF75EwghC0WSZs8q5UWWWv88jjdCridUEbtqfv5bsR1ziNiWMiZYIhaLJKAUZTdCJHXKNxO7ISIR1+RRi+8cc5CihAnIWePnPGH6vVH0+iLytYOmzdX+9x2ZfJ2OfFekHMJTnMKtTx11TTa1sJLeX9S6gIChIWc8YfKVWZUrsIUR9xe9PwIk7HIx64ICADK5fifX7uxW3vmSOLY0yVkWMgZf6hsQqSqsVNccXvR8ypM2az9koHWTJ5stvqzcT2/cWlclhTs/JagH+yRpwxVE5Bx8HSdCKqXt1+P3MtkX5zPr1uSNtmtGPBkJxMqjfzBxXXitJHsC6fjd3teTp4Uoqen/qSml06OjZ7fJhfPOMFCzoRHo2ln9T4XlaB4FUCVKZn10hC9XOQavSAlLXc/xbCQJ40kR0F+C0GcIlmVguLl/Hq1JIK4s1Bli1gj/Z6e2sce17ukJoWFPEkkPQoKwotVKShez6/Xfcf5+L0eexp89bigIDhjIU8SSY+C6o2/kS90rQKRoL14VcLf01N/bE4Eveam07En/bsYFxT9/ljIk0TSo6BaX9pGv9AqS7YbOb9eLj4nTwqxbl319lta/N1VqbDbvB47t2BQg6ILIgt5kkhDFOT0B+knc0LVUm5hnF+njJOof4eNHDs3RfOPouCMhTxJ8BfaHqugOFkt9bYVxvmN611VnL9baQhgnAg4IufKzjji1Fs7Desgqizp7+lpbFthnF+Vx1kPL/1p4vzdUt23J06orIC2w07dg35wRN7EqMwxb2mp9qKt24rKcw0r8vWyn7j7z2mOyIXgrBUmZais+uzpcd5WkGLq5hjCEE634hdnS0WShDFGTKBCDuCbAN4F8FM372chZzzTiOccVIQXJ8Fxe16SEu3G/a4hYpyEXMniy0S0C8AtAP9DCPHReu/nxZcZzzSykG9QixXHaVFht2NJ2sLNjC2BLr4shHgRwC9VbItpQtxM1jUyWRTUhGOtSTk/C2M08lm35yXMydewCHMRkrgveGIXpjfyANAPtlYYrzQ6WdfToz3qedRBWCC1/PpG9+dnrG79+rjYQSoI83hidO4Q9GRnPSEH8BiASQCTfX19oRw0kwAaLVAJOlOj1uec9u+lCMi6/TAKiNLkP4fp+cdofiFyITc+OCJndBopGc9mg/3DcnOhsBNFt8dit32nR9QFRHElzIKrGBV3OQk5FwQx9oTlCXrxbuWanuWy/WdUFY64WXLOblFht8fiZXHlJHvYQeLmXKv6DidhfsFO3b0+APwNgKsAVgC8BeA/1Ho/R+QxJ67+Y62yfJUReaMRmNtjqdXJ0frwsgKQG9Jir9Q71yq/w83kkXt5sJDHnLA9QbfiUksAVf5hOR1/Nuu8D3kM8n3yfDVyQQrSI4+JICmh1vdG9Xc4JhdAFnLGPV4i0jC/4I0IbCPU8rDdLpZcSyD9eOR+zndYF+gwvxNO+4qRr60SFnJGw80fmdty+McfD79XdZiWj9tJ1UYzb5y277QNv8dfT9xUCHBcbLkYZZqohIWccf9HZve+deu0JlVWAQj7jyXMaM9tVNdo9FfPK7f+bvyIU70LkyoBDlJAvaRsps1GWoWFnPH2R+b2jyaFt686bs+Xm/fZXYBqeeV2FymVk7BBRbBBWRqN2FEx8bVVwkLO+Psj85JpEZfbV+sf8uOPe/vD9nMH4yaDwqs11ajYuplbUCXAQUXkUU0QxwwWcsbfH5nTZxtZMzMM3ERwbvuguxH/RjMonD5n93qjdoEbka4l9l4i2qAsDbeBRFy+fwHBQs747+fhFFWqvH1VdTvsNoILI3pTueBxI+fHrfWj4sInt6Xa0nA7AZ9iEReChZyR+PkjC9pzVBnNuY3gwvDzvd4JeX1/vd+LF4tIbifoNgheSenkpVdSIeTj49qDSSkq/dU4ReReRchrHr9XkXZzEY5jHnYKJy+9wkLOxB+V4nHypPf0viDxIkJeLmhhTy6meCIxCTgJeSKaZj30kPY4fVp7yOeMd5SeO9WNtVQ2J9q9W5MeJ8JePX73bm2xh74+rbnXgQPO58vLIhpBrTwf9KrvjFISIeTNTiwvXLIT4eXLmmBevqw99yPmqsVj82bn12XHwrBwc77khfEP/gD4wAeAnh5tKbZaF516F79GL7a7d2v73Ly5/hiY6LEL04N+sLXiDRXHLbch75B9bzOoW2+VPmicJshqZV34GWu9DJe4HD+jBLBHnjxUiq9yIY/jZJgdcZkgq+XX+62sdDpG9rlTh5OQJ8paeeEF7cF4R5678XHt4ftcJqHZPlC9AAQQzSK6tc7LgQPBeN1B+edJIu6LJqvCTt2DfnAeuTdU3oko21bSbttPnrTvFxPWmE+edI7IZSTdSPTchB0AXZO076gLkAZrpVmJraUUF9uiHvWqFsMStiC69dUr/0+ZkHkihRcyFnKmealXHBSWr++muZbXC2MYPcaTSlLmcTzgJOS5iJ0dhgmeep5wWL6+TN2Tnnhfn5ZaKV/fvdt7el9fn5bKaCWT0R7WfTQTTucmbvM4CkjUZCfDNEStP9ywi1ysk69+BdYu9x4AymUt/lSR3w94nzSMwyRjMxU12YXpQT/YWmFCxckj7+lJh9UQdLMrr157nLz5lFlLcLBWSPtZuIyNjYnJycnQ98s0KadOAXv3Ateva897eoDjx9NpN2QymnRaIdLuAhqhv9/eopAVsn7fz7iGiKaEEGPW19laYdKNLI2XIg4At29HN56gCSK/30s++qlT9iJeazuMb1jImXRz4ACwvGx+bXlZez2NBOELu704yIum1+0wvmEhZ9KNUxR4+XI6K/2CaHbl9uJgd9Gs9X5GGeyRM+nGya+VtLZyVz83nDrlnDYpcfLnAeDkST7HCnDyyFnImXQjb/edIkWAJ+FUwZOcgcOTnUxzYrQanOBJODU0U952zGAhZ+KP3+ISWYTjJOY8CacGXowiMljImXijciUijhiDR3XlalDEofJUISzkTLxRmT4YdMSYMnFILUEsUxgxPNnJxJsgKhX9YpfBAVRPqnJGTDxJ8KQsT3YyySRuKxE5RXN798az8IjvEqpJ4cpJSoSciB4mol8Q0UUi+lMV22QYAPHztZ2sHmMLACNRikMKLQQlxC04UIBvISeiLIC/APBpAPcB+D0ius/vdhkGQPwyIbwKc5Ti0GztCdwSt+BAASoi8u0ALgohZoUQdwB8C8CjCrbLMBpxyoRwEuaenviJQwotBCXELThQgAoh7wXwpuH5W6uvmSCix4hokogm33vvPQW7ZZgIcIrmjh+Pnzik0EJQRpyCAwWENtkphHhWCDEmhBjbsGFDWLtlGLXUiubiJg4ptBAYe1Ss2VkA8GHD8w+tvsYw6aSRtTWjoN4aoUxqUCHkrwLYQkQD0AT8dwH8voLtMgzjlwgvOkIIEJHjc0Ydvq0VIUQJwJcA/COA1wH8LyHEz/xul2GY+GMtKJTP5+YO4eLFJ/TnQghcvPgE5uYOhT3EpkBFRA4hxPMAnlexLYZhksHc3CGUSgsYHj4KItLFOpfrQqm0iELhOABgePgoLl58AoXCcfT27uXIPAC4spNJP1zd6Btr5F2pVFAqLaBQOK5H3lKsS6VFDA0dQW/vXhQKx3H6dEYXcSn6jFq418oqDz2k/fvCC1GOglGO3cIS3APFE06RdzbbhXJ5LfIGYBJrIQROn16LFcfHKyziPuFeK4p46KE10Y87SRprYMS1ujEhdwlCCMfIu1zWIm8jVrE3YvTMGbUo8ciTjBS606fNzzkyTwlxrG603iXIHihA7O4SiAjDw0cBAIXCcT367u3di6GhI5iZ2W96/8WLT+ivG+0UKf4A2F4JgKYXcrckSfCTNNbA6euzb1ka1x4oMRNyYE3MjRZKPbHOZrtMNou8GORy3SziAdD0Qi7FranFLs089ZS9R75a3RjJ772RuwQ3q9gHhJ1NMjOzH7mcs1gPDBwyZafIn1tFnHPN1dD0Qu6WJAl+ksYaOHGsbvRyl3DqlNbr3NgmN0QrxuiJWyNvaa84ibVVkK3PndMXtQsB4x4W8lWaWuzSjk11Y6T2U527BB27jBtJSFYMESGX63aMvDOZTNX73WCcRAWCyzVvloif0w+ZpsQq5OPj2r+hXdDdWCVOS5JJQlzuLghBNEb7EpW55mmM+J3SDzkiZ5qSyO0nq+Uj0yGNYl4vsybECVsnm8SPwNtNohpF14+Ye4n40xC1cx45w0SBm2XYagl1DNrR+u2n4pRrXqlUfPdlkReJetWlaekJw0LONDUvvBDR/IibQiW7fuKAthrRs89C/L65yWiYNmmtQqFSaaHuWIzvb2sbAQC0tY2gUDiOqaltq9u5YRJYr8dn9PMlxojfeAwXLuzzfAxxgq0VJpZYK1JTNxntJgWxRsbN3NwhlC4+EZn/S0QYGjoCIYSpUCifn3DlcRMRbt2aRlvbCEZHJzE7+6S+jWJxGi0tG00ifuHCPiwtncUHP/iw6+NzivhlU6/h4aMYHj4KIQTm509gfv4EALU+fVhwRM4wUeB2GTabVYf8RsMqmJs7hPPnH4RV65aWzuLSpcN1Py+EQHv7CIrFaczOPllV6r+y8jbm55/BxYv7cOHCPszPn8DS0llTlF5v+0ZPfHy8otss1659Tz93AKqOIWkiDnBEzsQMazaJ9fXUROZuUxBtqFU2bxShoCbxhBBYWbmBpaWzWFo6a/rZ0tJZdHbuqLsvp2OQtLVpIl8onNBf6+2dwPDwMVfHUCttMpvtQnf3uO1+AS1qT5qYc0TOuIIbcCnG50rutfxfwPskntMCEXY/k/tqbd3qsC13fr20Z4zs2lVGb+9eFIvTVe8fHDyi20hO4zIyMHCoqkBpePgoBgcPV527fH7CFLUnrcEXR+RMrLCmBVpfTxU+lmFz8n+lQHkptpmbO6h7xkSESqWyWoLfrW/Lzotfv/5RXLnymmkM+fwEcrlu/fP9/Qcd7wqEEFVNt2Zm9mNw8Gm8++63sLLyjulnL7+8Effe+x8B3EI224ly+abtnIC1NYDdc+u5k0NMak8YFnKmJtyAK37UKpsHoE/iAbWtFwCYnT2I69e/p0fAQ0NHMDW1DcXiNHp7JyAE9ElA8wVhApWKfcRaKi1gfv4EOjp2oFS6gaGho8hkMqY+5gMDh6pK/WUTroWFF3QRJ7obQry/ut1rKBT+G4Aycrn1KJWuVV2Y6hUB2Z+7fasWDunnLkkiDrCQMzGFLxTO1C6b79JFyKnYRiKEQLm8iGJxWk/9k+9vaxvB0JD2fiKyXBA0Eb969Rn9+Y0bp7G8/Jou+ps27UEmQygUTmBh4UVs2zalC3Vb2wiI4HgMt25No6NjO27fnkWpdM1y9GUAmqi3tY1gcPBpU4S/snLD4cKzJvjV+z0GgBIXhRvhEv0EE2Z0zJF4/LBaJLOzB1EuL+qCuBZpathF5HZl8oDmVcs+KnYr/Vy6dBi//OU/oLNzB4aHj0EIgRdfzJreI4TQo3uJnMQ0Ntyys16EEKhUKnjpJedYs6/vP6NcvmmJvvfh5k3zJKzTcSexmpNXCGKYlGEXXcviFqOI5/MT6O2dsJ3Es5twBDSvWhbN2Hnx/f0HMTr68mo0C9sFJogI27ZNmV43ivjMzH5cunTI8fis+7Vy/frf2aRgnkBn5w7T+xqxSmpN/sYRtlYSSBS+NUfi8caazifJ5yewZcsx+a4q+6BSqWBqaptpW9JmAYTukTt58QAc/Xotyq0e6+Dg0yabZfPmgwYPfR9yuW5cvfoN3LlTQDbbg3L5etU2crn1tpaQ9PWNWNMJ6/noSWy2xRE5kyg4DdIZu5TELVuO6fbF8PBRkxDJrBEpiMbUv7a2EeRy3WhpuafKx+7t3atfEJz8+nx+AktLZ1EonNBL8CUvvZRDoXAcLS0bUSxOY2bmCV3EC4UTuH79HyBECQBQLl9HNttj+nw226N75D09j5h+Zrzw2KUTasVUN/TXtL4u+/RiqnK5HHmxVSOwR55gHnoImJ4GRkaaJ2Jmr96ZRtrCyuhzaOhIVWbJ4OBhfbv1/GS798zNHdIzYnp792Jw8GmT553LbUKpdLVqTC0teezceQVnznwIKytv66+3tm5FV9cutLR0o1RaRDbbhUrlpul429u3o7NzJ4aHtUwZYyqljLZXVm6ACKb5g46OHbj//h/rKx+VSouBtdf1A7exZRINp0HWxk1Kop0IuVmSrd5KP07vGRw8DCKgu3vcdqFmYMX2WNat+zUQER54oGCaQB0bOwciwqVL2nbL5UXMz59Ae/t2CHEHKyvv4tatV9DZuUO3QsrlReRyXfpxytTITZv2mPbZ3r4dMzOaxy49/FoZP3GDhTyhSCFbXNTELe7CFvfxJZ16K/nUK5ev9dyYSWL3by0GBg7rUbH0xKV1o2WzZCFTCrWxrsfo6CSEEDh3zhx4Tk1tw+jo5Kr1oeWp5/N7sLj4EorFaeRy6wEA7733v7Gy8rYpQ0aO1dgky8haKuVe24tO3Mv2WciZwFAp3pEvBJEA3C547AVpvWSzXSiXF3WRs1ZW1rJfMpkMslltoeaBgT/H3NwfG6Ldsml/pdI1nD//AN5/f073wUdHJ/Hyy5qffu7cGEZHNVu2UDiupxmuXRig2zFSxKvvMJyP1zgR6+XOJmpYyBNKUoRtelobI1si4eDGBnGLscuiFMqFhdN69CsrK+18aKvvXi4v4r33nsPCwmmMjk7aNquSvP/+ZZRK1/TofHb2Sf15NtuJbDZbVey0bduUyYaR2AnvzZtnq94nmZ19Ur/oeL2ziRIWckY5Vj97elrdtvkCEC7WlEYZ9cpo2VhaL0Vdiv/CwmlTRSfR3bhz5wpefnmjZS+ETZu+iMXFl7C8rPVuaWvbimLxNX1ytLd3AoODR5DNZlEulzE7a7Y+rCmUEqMlIucRlpbOmipPAWMl6nFTsRKg5s4maFjIE04ShG1kpPmya5KOMZfaGv1KisVpPQo2RrBDQ0f0yF3+XEb0RHdbyu4zACq4evVraG/fjtbWrSiVrqFYNDfjqlQEMpkMZmcPYn7+66sR+r3YufMtnDmzafUCcxeAX5l8eKslYpxHuHTpsJ533tJyD/r7D0Lm2suqVkmcRRxgIWcCwM724dzv5GC0VJyKeqwYI9ZMJlNldRSL08jnJ0DUgUJhref6gw/exvT0diwvv4Zbt17RXzd63oA2GSmEwNLSj/QLQan0Ds6dG1t9ngHwK7S03Iv7738V586NIZvtNOW8A1gVa02YBwYOoVKpmNoExD3ydoKFnAkFjsTDw28fkVrZHVJgpUcusbMw7Ma1sPD3ptd+/OO7AGg54tJWAaB3X6xUBK5d07JQ3n77a1XbXPtMBQCwYcPvYGZmP5aXX0NHxw7dpwfMdxlyPEZvXx57EvFV2UlEnyOinxFRhYiqktSZ5KGycjKyhY2bGLcLSjj1EpH/EpGhtF+jt3cCo6OTaGsbQal0DS0tG9HX9xVT9WS5XMa5czsxP3+iauGJq1efMYm1EevrbW0jeofF9es/5/r45+dP6KmEnZ07TJ0Rg6jYjEtPFr8R+U8BfBbAXyoYC8MwPjCKFeDcxtWpl8itW9Nobx/RJ/ouXtxn2v6NG6dXbYh1AICenn+rWxUAsLR0HufPP4j3378EQBPn1tatuHPnbZRK5kUirBE4oPWFATQxLhanUSxO65OQXuntndDb8ALul8fzQpx6svgSciHE60Byb0eYNbhyMvm4EataYi8nBxcWTqOr65OYn9ci29bWrVhZeRvLy6/pvnd7+3YUi6/g3LkH0NW1EwMDT+P8+TF9kpKoDUIUdbHOZtejXF6zYpwi8+HhoyY7x5hZ0tb2cQArtsvAWZGi2tJyjynPvV6Pdre4vWiGBXvkTFPQLBememJlbFtrFfvBwadx7tyYHg0Dxv7hE6beJESkF+PcuvUKFhZeNGWaCFE0jatcvobe3glks92Yn//Lqgg9l7sXuVxXlbe+uHgGmUwbKpUiOjt3Ipfrxu3bV1Cp/BKAljb47rv/s6pDorwI5fMTep92p4pNo4+ujb2+CAcR4fuhrkdORD8kop/aPB71siMieoyIJolo8r333mt8xEwgSD97fFx7pN3fTmsXRaf+4dK7nZs7hJmZ/VU9yLPZLmSzWYf+4dWtYW/fnqt6Xz0qlQpWVm7oIt7Scq/e86RUegfXr/+d3jtFdky8desVVCraReH69e+gXF7Exo1r65xevfqMLuLZ7HrT/lpaNmJo6Ijep31qapt+5wGsteudmtqGSkWbLK23SLURu26TUWW91BVyIcSnhBAftXl818uOhBDPCiHGhBBjGzZsaHzEDOMBKdinT2uP7m7tkUasjbOsbVyNxTqTk6Omz1679l2Uy2WbxlbQs1dkq9t8fgKl0ruuxyUnPa9e/Zo+EQkAKyvv4OrVZ/SFL3K5LuTzE+jq2ml7Ybhzp4D5+ROYn39G99MlmzZ9CQ8++LbptXz+Mb0KVLbnBdZK9+XErdZKd7/nCdB6F80wUdLGloheAPDHQghXvWm5jS0DhGN3WL3/7Gpqc3m1xcf4ePBjCJN6E3CVSgUvv5w3rVAvJx5zuV9DqfRuVQ53Lncv1q3biOXl1/Sqx1dfvR+3b//EcRytrVvR3b3L5LM7+eK7dpX1XG65xJtdC1unzwNaJSjRXaZcdOvcgHW5OiIyNfSy+5wTtbpNBmmvBLLUGxH9NhG9BeABAP+HiP7Rz/aYdODXtggiBbKrS3uUy2sinkYGBg5VeeLGBSWICHff3W/6zD33aFezSqWIXO5evSryk58sraYaaqLf0bEDhcJxvPhitqaIA7J17VG0tY0gk2nHnTtvO753amqbKf3x3Lkxk4gDzpOjcozFolZQtGnTHts7EafIOZPJNGSPOHWbtBYghYXfrJXnADynaCxMxIQ1IRiHDJmurvD3qRqnwp96jbPa27ebFieWq/hoAv4vUS5rnQZnZ5/E6OikPgEqV/2pR0vLvcjlujE7+6SheOgdtLWNoLPzk7q90tq6FUIIFIvTeota82euVW13w4Z/p0f5APReLC0tG3HXXX34yEeOWxpdddXsZuinZW0Q3SYbhbNWUkoU4uhXoL1+3sv2FxbMn0k6jeQwy/dcvfoMNm3aY/Kr5aTm8PAxVCoVzM4+acrG0Mrrq7eZz+/B0NBRzM7uR6FwAplMK1ZW3sGVK/8FAFYj8g7cvv26KRumpeVeZLN3YevW/4vz5z+OYnFaF2V5UcnnJ0ypiC0tG03+uPFnO3e+hUwmYyuqc3OHbLsZZrO1Rd5tZF6LsNIQec3OlODHjrBOCD70kDYhGJToxSFDJslZOY1WKUo7IJ+fsC2yGRrSBG529knkcl2Wz2qRe0vLRuTze7BuXR6AViR0/vyDIOpELrcelcqy6XPF4jQqlaWq6HrDht/B0tIrOHv2w/jgB3/L9LOent+yvXBobXMnqqwQOWbrOZICOjBwyLab4eDgYaX2iNuq2iDgiDxlRGlb+O2R7vbzfo4xqeJtxE8Oc3//QX2RY+uk5tTUNt36sC6YfPPmGX3CcX7+GbS03AsAuldeqfwKpdI1ZDKttmJutUoWFrRf3srK23jzzT8zvb9UWtQvHGuRsjbmSkVgZuYJfYFlYxtdAHpXQ+vdilNPFVX2SNQFQizkCccqal1djm91xCigsnd4GEvIpUFUo6LRKkUiQjbbbSr0kT1NpPXR2rpV/3f9+kf1POx8fg+6unbh6tVnTFkvwNpkpFXEJVr/8q160ZDd5OWuXWVdlDs6dujRt1bEdBSyxezCwgt6L/RMJqO3zb11a3p1X7UFVZ6HWufIK1EXCLGQp4yREfPzRsTS70IQfgW63ueTsjpSkDjlMLsRDeuiyEIIk1/e0/MIAE1sb9y4G93d/xptbSOrPbsPmd5rpL394xDCvoS+rW0EHR2fqOozbsRYqJTLdaO//6CpN8xaD5gbKBRO4Pz5BzE6+jIuXnxCzw0H4Fi5Ojx8FHNzh1AuLwbSH0VlCwCvsJAnHCdRa8TftkbmvBBEPKmVwwy4E4+BgcOYnT2I8+cfRLn8vulnRqujs3MHyuVFFIvT6OrahQsX9tYY14ppYWWjdWOc6DQi1+SUk6sATCXzdpaFnAJYWjqr54bn82veuWajmG9NZVve69e/p49Dtf3h5+LqFxbylNKIAFttGrneZlzFPK7jChqnHGbA3bqS0l4olxf0dMLW1q0ggilizuf3YHj4mP7cmCWSzfZU9TeR4t3T84ge7c/MPGHq0QJonQkXF8/g1q1XVtMc9+sTrdbVeZwsC2tGizxkYwMwIxcu7Fs9vrWVg1TaHyourn5gIU8JQYia1aZh4oNxks7Y2c/43A7zEm7HcPPmWSwtnXUouFmbANyy5ZhJOFtbh7G0dB2trVuRydylN9EiWmeyKKwJNJq3rY3zjTf24tatV5DL3aMX5tiN286ysL6tUDhh6JK4lr4oJ03l2GV7W+PqRSpE1u/F1S8s5IwOe8/JolZvcTvP196mcG7RIcVvePgozp9/sGpbbW1b0dPzCPr7D+n9y3O5e0zjkNklw8NHceHCPszPn8DMzH4MDx/FRz5yXD8O479y+8ZFIS5c2Gfav5bRouW9y+1K5B3BWiS/9rOhoaMNFwDVI8oCIRZyhkkoXlPenGwKO2RUu7R0Fhcu7NMtmI6O7ahU7ug9Tbq6dmFm5gnkcvdgePiYScSsEeqWLcf012uJm3VyUwp1e/t2bNt2BufOPYClpbMQQgq++fOyZS2AKs96amqbPjEahP1Rr6o2KFjImSo4Ek8GjaS82dkUdnR0fALd3buQy90DQK7eI0zl8e3t20FEer63lUYiVOPFaWHhNEZHJ/WLiBB3MDt7EKOjL+PChX3I5bpX0xXXovM1YRYQAqY7AqN/Li8SYdofQaKk+6FXuPshw6jDqbOf03uNUSgAfOADHzM1wWpv346enk/r6X/ycwBM+5HIwhzjJKUfKpWKHjlLjJkw69ev2Tm1rCUADj/rwsDAYdM5SYqIO3U/5IicYRKMl5Q3u8yKN97YW5UX3tW10yTixu3aIasr6+ViWydnrf9KMpkM7r//VfzoRy36a8Z0xu7ucf31WlG/sfhH/t9uQjgpIl4LFnKGSSheU96svjUAvedKR8cOvbhG+zxVRbKFwnFs2rQHN2++VFXaL31np+hWRs65XBdKpbVl17LZLpTLi6aLwOzsQVy//r2qbRi9betxGc+JtYhIK9m/ASGwWtR0MLJFkoOChZzxBGe0xIdGUt6sEWwud4/uLzt93thsy5iLbSz2MfrOVoy+t/zcwsJp03bkRcBYtGPXynZg4M/rplYODR0x+eyZzDosLWmTs/n8hN63JYpFkoOChZxhEkwjE4rGn7n9vHzfpUuHdU/cmIvd0/OIo0dunZQFYLoIGC9ERISenkfwq1+9VSXiAHD+/MexbduU/l6JNYNH9l+xVpOu5ZNHs0hyULCQM66wVn1yZB4f/Ka8efm8XC7OmotdLi/WjG5rZctYBXVw8DCEqODNN/+r/lo+vwcLCy+iWJzGxYtPgEi7mzB2M3SbWmm3z6TD/cgZhqmL7LVtXOOyrW0EfX1fMS2r5pQFZzcpK7F+Tlu3c8n0nvn5Z9DdPa6vUlQonKjqvW4U83q4WSTZ+vMoMvzcwhE54wqu+mxerLZFLtdlyiDRim+EyVe3VmYac7itjbWMk7MATBO4Q0NH9ApRY/WmnTVid7FoaxtBV9cn9fz3tbL92kVAjazAFCUs5EzgsPjHG6e1PyV2HjewJqaXLh2GENqiFfLzRtHTJku7Vlfe6UJ39zgGB5/G7OyTyGY70d09rl8EhBC2FaGAuWmXk4jLC0A226VPmmYydyGf3wOA9KwV2dvc7rijXiSiEVjIGU+wGKcLt5Gnncctxb1UWsD8/An9PVbR09L/FvWiIWnPLC2dR0fH/frrxn1Xpxiax23NlbfL4OnvP6jnt8uLjHy/dfvW8yAnS6NYJKIRWMiZwOAJ0njjJfKsVXhUq00AsLZijxACW7Yc0z321tatjvs2jlF7vbrUXn5OjtGagVOvo2Kt8zAzs78q4yWuIg6wkDNM0+K2V4ubwqNaK+MMDx/FzZtnTD639Kqz2btr7ttrrnwjGTxO58G6pmlYi0Q0Ags5Exg8QRp/nCyTWhWhVjEFqsv3raLX2blDb36lbXOtFa3x9Vo57F5y5b1idx6C7JKoGhZyhmli3PZqcRJT+f5a0bq2H/N+pYhbX3eKev3mytfDKeMlKV0SOY+cCZwXXuBoPI5YLZPx8UrNnHA7MXWK1rUMlW4A0NMHtVa4a9y4cVpvM1tv30HidB6KxWnMzOw3rcAUx9RDgCNyhmlaGunVYkc968PYp8VIqXRNXzQ5yqjX7XmIYyQu4X7kDNPk1MsjV7F9Y6MqzX7ZB+OCEHaFRGET9HlQAXE/coZh7Ajaf9Yi3nssEe8x2BXlRCmcQZ+HIGEhZxgmcKJcmLgZ4MlOhmFCIckRb9xhIWcYhkk4LOQMwzAJx5eQE9FXiejnRPQTInqOiLoVjYthGIUkqbc24x2/EfkPAHxUCPExAG8A+LL/ITEMoxK5KIQUb1kAMzd3KNqBMcrwJeRCiO8LIUqrT88A+JD/ITEMowpjZz8p5rKK0brCDpNcVKYf/hGAv3X6IRE9BuAxAOjr61O4W4ZhnHDb4ZBJNnUrO4nohwA22vzogBDiu6vvOQBgDMBnhYtLPFd2Mky4CCFw+vTaDfj4eIVFPIE0XNkphPhUnQ3/IYDfBPDrbkScYZhwcdvhkEkufrNWHgbwJwAeEUIsqxkSwzCq8NrhkEkmfj3yrwG4C8APVq/sZ4QQX/A9KoZhlKCqwyETb3wJuRBiWNVAGIYJBu5zkn64spNhmgDuc5JuWMgZhmESDgs5wzBMwmEhZxiGSTgs5AzDMAknkjU7ieg9AJdXn64HcC30QcQHPv7mPf5mPnaAj7+R498shNhgfTESITcNgGjSruS0WeDjb97jb+ZjB/j4VR4/WysMwzAJh4WcYRgm4cRByJ+NegARw8ffvDTzsQN8/MqOP3KPnGEYhvFHHCJyhmEYxgcs5AzDMAknFkJORH9GRD8homki+j4R5aMeU5gQ0VeJ6Oer5+A5IuqOekxhQUSfI6KfEVGFiJomFY2IHiaiXxDRRSL606jHEyZE9E0iepeIfhr1WMKGiD5MRP9ERP+8+r3fq2K7sRByAF8VQnxMCDEC4O8BfCXi8YTNDwB8VAjxMQBvAPhyxOMJk58C+CyAF6MeSFgQURbAXwD4NID7APweEd0X7ahC5b8DeDjqQURECcCTQoj7AOwE8EUVv/tYCLkQ4qbhaRuAppqBFUJ8XwhRWn16BsCHohxPmAghXhdC/CLqcYTMdgAXhRCzQog7AL4F4NGIxxQaQogXAfwy6nFEgRDiqhDi3Or/lwC8DqDX73b9rhCkDCJ6CsC/B7AI4F9FPJwo+SMAfxv1IJhA6QXwpuH5WwB2RDQWJiKIqB/A/QDO+t1WaEJORD8EsNHmRweEEN8VQhwAcICIvgzgSwAOhjW2MKh3/KvvOQDt1utUmGMLGjfHzjDNBBG1A/g2gH0WR6IhQhNyIcSnXL71FIDnkTIhr3f8RPSHAH4TwK+LlCX3e/jdNwsFAB82PP/Q6mtME0BELdBE/JQQ4jsqthkLj5yIthiePgrg51GNJQqI6GEAfwLgESHEctTjYQLnVQBbiGiAiNYB+F0A34t4TEwIkLbG3jcAvC6EOKJsu3EI/ojo2wD+BYAKtPa2XxBCNE2EQkQXAdwF4PrqS2eEEF+IcEihQUS/DeAZABsALACYFkL8m0gHFQJE9BkAxwBkAXxTCPFUtCMKDyL6GwAPQWvj+g6Ag0KIb0Q6qJAgok8AeAnA/4OmdwDwn4QQz/vabhyEnGEYhmmcWFgrDMMwTOOwkDMMwyQcFnKGYZiEw0LOMAyTcFjIGYZhEg4LOcMwTMJhIWcYhkk4/x+M7bm+pC+lugAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "draw(X, Y2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "20574922-4111-43ea-aadb-99bbe5cd640e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "dc136240-a72e-4025-b345-5075b34a8f54",
   "metadata": {},
   "outputs": [],
   "source": [
    "datapath = 'Rice_Osmancik_Cammeo_Dataset.csv'\n",
    "data = pd.read_csv(datapath)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "7223ff2a-084d-45fd-97a5-b455e939a006",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(7, 3810)\n"
     ]
    }
   ],
   "source": [
    "#处理数据，将数据存入矩阵中\n",
    "attributes = data.columns\n",
    "Y = list(data[attributes[-1]])\n",
    "for i in range(len(Y)):\n",
    "     Y[i] = 1 if Y[i] == 'Cammeo' else 0\n",
    "Y = np.array(Y)\n",
    "Y = Y.reshape((len(Y), 1))\n",
    "attributes = attributes[:-1]\n",
    "X_uci = np.zeros((len(Y), 0))\n",
    "for attribute in attributes:\n",
    "    X_uci = np.hstack([X_uci, np.array(data[attribute]).reshape(len(Y), 1)])\n",
    "X_uci = (X_uci - np.min(X_uci, axis=0)) / (np.max(X_uci, axis=0) - np.min(X_uci, axis=0))\n",
    "X_uci = X_uci.T\n",
    "print(X_uci.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "a034967d-8933-4439-8edf-e5af09b65115",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████████████████████████████████████████████████████████████████████████████| 50/50 [01:16<00:00,  1.52s/it]\n"
     ]
    }
   ],
   "source": [
    "k_uci = 2\n",
    "alpha_uci, miu_uci, sigma_uci = GMM(X_uci, k_uci, 50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "515dc885-3125-4f81-b64a-43d74368eee0",
   "metadata": {},
   "outputs": [],
   "source": [
    "Y_uci = np.argmax(GMME(X_uci, k_uci, alpha_uci, miu_uci, sigma_uci), axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "806a9775-5ee7-4e23-87b4-027c6a8a8752",
   "metadata": {},
   "outputs": [],
   "source": [
    "neg = X_uci.T[(Y==0).flatten()].T\n",
    "pos = X_uci.T[(Y==1).flatten()].T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "2fff6178-e6bd-477f-b8df-584ae7829b35",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[array([0.5819303 , 0.67780122, 0.64231117, 0.60893771, 0.72501644,\n",
      "       0.59523796, 0.42351801]), array([0.35194363, 0.37135988, 0.33092909, 0.51961286, 0.57993619,\n",
      "       0.35834966, 0.47405267])]\n"
     ]
    }
   ],
   "source": [
    "print([np.average(pos, axis=1), np.average(neg, axis=1)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "8eadc06d-58a8-4e51-bdb2-809df9095a5c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[[0.5475837 ]\n",
      "  [0.64250895]\n",
      "  [0.61547313]\n",
      "  [0.57793027]\n",
      "  [0.73206365]\n",
      "  [0.55975932]\n",
      "  [0.42061587]]\n",
      "\n",
      " [[0.35076753]\n",
      "  [0.35907053]\n",
      "  [0.30920285]\n",
      "  [0.53724553]\n",
      "  [0.54979478]\n",
      "  [0.35724204]\n",
      "  [0.48500971]]]\n"
     ]
    }
   ],
   "source": [
    "print(miu_uci)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "a2644c36-7680-4d68-93af-2fc716258f03",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.01282604  0.01236099  0.00981268  0.01075683 -0.00159672  0.01300398\n",
      "   0.00157676]\n",
      " [ 0.01236099  0.01372367  0.01201759  0.00833073  0.00124411  0.01265615\n",
      "  -0.00093275]\n",
      " [ 0.00981268  0.01201759  0.01215119  0.00432064  0.00396531  0.00996339\n",
      "  -0.00097965]\n",
      " [ 0.01075683  0.00833073  0.00432064  0.01241889 -0.00579671  0.01093417\n",
      "   0.00296345]\n",
      " [-0.00159672  0.00124411  0.00396531 -0.00579671  0.00613934 -0.00162518\n",
      "  -0.00250086]\n",
      " [ 0.01300398  0.01265615  0.00996339  0.01093417 -0.00162518  0.01324884\n",
      "   0.00145011]\n",
      " [ 0.00157676 -0.00093275 -0.00097965  0.00296345 -0.00250086  0.00145011\n",
      "   0.05109478]]\n",
      "\n",
      "[[ 0.0084091   0.00901939  0.00680862  0.00822557 -0.00164826  0.00854754\n",
      "   0.00033203]\n",
      " [ 0.00901939  0.01132993  0.00962202  0.00675856  0.00185804  0.00925804\n",
      "  -0.00181338]\n",
      " [ 0.00680862  0.00962202  0.00997399  0.0024199   0.00589514  0.00689157\n",
      "  -0.00217744]\n",
      " [ 0.00822557  0.00675856  0.0024199   0.01219901 -0.00863117  0.00842571\n",
      "   0.00252458]\n",
      " [-0.00164826  0.00185804  0.00589514 -0.00863117  0.01237718 -0.00175478\n",
      "  -0.0039804 ]\n",
      " [ 0.00854754  0.00925804  0.00689157  0.00842571 -0.00175478  0.00872826\n",
      "   0.00022196]\n",
      " [ 0.00033203 -0.00181338 -0.00217744  0.00252458 -0.0039804   0.00022196\n",
      "   0.03957484]]\n"
     ]
    }
   ],
   "source": [
    "print(np.cov(pos))\n",
    "print()\n",
    "print(np.cov(neg))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "d8b899ff-b885-422e-8eec-f8eeaf7457eb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[[ 0.01798309  0.01805533  0.01436265  0.01516998 -0.00214361\n",
      "    0.01835414  0.00177098]\n",
      "  [ 0.01805533  0.01970937  0.01643017  0.01371333 -0.00010292\n",
      "    0.01856315 -0.00034671]\n",
      "  [ 0.01436265  0.01643017  0.01506527  0.00914278  0.002132\n",
      "    0.01467177 -0.00029386]\n",
      "  [ 0.01516998  0.01371333  0.00914278  0.0153959  -0.00509774\n",
      "    0.01551872  0.00272189]\n",
      "  [-0.00214361 -0.00010292  0.002132   -0.00509774  0.00448189\n",
      "   -0.00220473 -0.0018953 ]\n",
      "  [ 0.01835414  0.01856315  0.01467177  0.01551872 -0.00220473\n",
      "    0.01879786  0.00166176]\n",
      "  [ 0.00177098 -0.00034671 -0.00029386  0.00272189 -0.0018953\n",
      "    0.00166176  0.05224083]]\n",
      "\n",
      " [[ 0.00902782  0.00934428  0.00665067  0.00942381 -0.00271244\n",
      "    0.00918086  0.00058912]\n",
      "  [ 0.00934428  0.01069073  0.00807124  0.00871834 -0.00077724\n",
      "    0.00958628 -0.00098971]\n",
      "  [ 0.00665067  0.00807124  0.00721632  0.00465475  0.00221192\n",
      "    0.00674195 -0.00088572]\n",
      "  [ 0.00942381  0.00871834  0.00465475  0.01221866 -0.00711824\n",
      "    0.00963994  0.00178865]\n",
      "  [-0.00271244 -0.00077724  0.00221192 -0.00711824  0.008671\n",
      "   -0.00282742 -0.0025253 ]\n",
      "  [ 0.00918086  0.00958628  0.00674195  0.00963994 -0.00282742\n",
      "    0.00937347  0.00047688]\n",
      "  [ 0.00058912 -0.00098971 -0.00088572  0.00178865 -0.0025253\n",
      "    0.00047688  0.0357003 ]]]\n"
     ]
    }
   ],
   "source": [
    "print(sigma_uci)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "id": "e8feef49-8ada-48bb-a266-42355971469f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.42782152230971127, 0.5721784776902887]\n"
     ]
    }
   ],
   "source": [
    "print([pos.shape[1] / len(Y), neg.shape[1] / len(Y)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "id": "7f7f7c93-73e2-4279-9523-556cef9c43d4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.5059002 0.4940998]\n"
     ]
    }
   ],
   "source": [
    "print(alpha_uci)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "393a1f5d-f10a-4819-8f23-24a74771fcf3",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
